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Abstract 


Efficient and accurate specialty finite elements methods to analyze textile compos- 
ites were developed and are described. Textile composites present unique challenges to 
the analyst because of the large, complex "microstructure". The geometry of the 
microstructure is difficult to model and it introduces unusual free surface effects. The 
size of the microstructure complicates the use of traditional homogenization methods. 
The methods developed constitute considerable progress in addressing the modeling 
difficulties. The details of the methods and attended results obtained therefrom, are 
described in the various chapters included in Part I of the report. Specific conclusions 
and computer codes generated are included in Part II of the report. 
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Chapter: 1 


Macro Finite Element for Analysis 
of Textile Composites 


ABSTRACT : The analysis of textile composites is complicated by the complex micro- 
structure. It is not practical to account for this microstructure directly using traditional 
finite elements. A new type of finite element was developed to efficiently account for 
microstructure within a single element. These new elements, which are referred to herein 
as macro elements, performed well in initial tests. 

INTRODUCTION 

T wo OF THE major obstacles to widespread use of laminated composites in 
high performance primary structures are the low strengths normal to the 
lamina and the labor intensive fabrication processes currently used. There has 
been considerable research aimed at developing tougher resin systems to enhance 
the through the thickness strength. Also, robotics are being developed to reduce 
the labor costs. Of course, there remains the question of whether laminated con- 
struction is the optimal form. 

Several alternatives which are receiving attention are weaving, braiding, stitch- 
ing, knitting, and combinations of these. These various forms are referred to as 
textile composites. Approximate analyses have been developed for predicting 
moduli, but these analyses are far too crude to predict details of the local stress 
field (l— 3J- Very little detailed three-dimensional analysis has been performed. 
These studies, which used 3-D finite elements (4-8], required tedious modeling, 
many simplifying assumptions about the material microstructure, and only con- 
sidered very simple loading. The computational challenge is obvious when one 
examines the schematic of a simple plain weave in Figure 1. The resin pockets are 
removed to show the fiber tow architecture. This tiny piece of material, which is 
only about .28 mm thick and about 1.4 mm wide, is in fact, a fairly complicated 
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structure. If four mats are stacked to obtain a thicker composite (still only about 
l.l mm thick), it is obvious that the number of elements required becomes intol- 
erable very quickly even for a coarse mesh (see Figure 2). A variationally consis- 
tent and organizationally (and computationally) tolerable procedure is needed for 
analyzing textile composites. 

The objective of this paper is to describe a displacement based finite element 
which accounts for the spatial variation of material properties within a single ele- 
ment. This is in contrast to the usual choices of either adding more elements to 
account for microstructure or using averaged material properties within each ele- 
ment. The performance of this element is very similar to that in Reference (9), 
but the formulation is totally different. The formulation of this new element will 
be discussed first. Then several configurations will be analyzed to evaulate the 
performance. For simplicity in the discussion, only two-dimensional configura- 
tions will be considered. However, the approach is general and can be extended 
easily to three dimensions. 


THEORY 

To simplify the discussion, a rectangular element with multiple layers of 
materials will be discussed first. Such an element might be used where the tows 
are straight or for ordinary laminated composites when there are too many 
lamina to model each individually. Then, microstructure of arbitrary shape will 
be considered. 

Consider the four node rectangular element in Figure 3 which contains three 
lamina of composite material. To facilitate the following discussion, the element 
will be referred to as a macro element and the subregions (lamina) will be refer- 
red to as subelements. The displacement field within the macro element is 
assumed to take the form 


u(x,y) = N { (x,y)u t 
v(x,y) = Ni(x,y)v; 


(D 


where N,(x,y) are interpolation functions and u : and v, are macro element nodal 
displacements. In Equation (1) and subsequent equations Cartesian index notation 
is used. In particular, a repeated subscript indicates summation. In Equation (1) 
the summation is for the range i = 1 to 4 since there are four interpolation func- 
tions for a four node element. The assumed displacement field is referred to 
herein as single field because a single approximation is used through the entire 
macro element. In contrast, a multi-field approximation would use approxima- 
tions which are defined within a single subelement. The stiffness matrix can be 
calculated using the familiar formula 
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(2) 


K.j = 


Si B ~ D ~ 


B nJ dxdy 


where B n , and D„„ are the strain-displacement and constitutive matrices, respec- 
tively. They arc defined by the following equations 


= 0 -/ 4 / 

o- = CL„c„ 


( 3 ) 


where q, = list of the element nodal displacements. 

The complication that we have is that the constitutive matrix D„ n is now a dis- 
continuous function of position. However, because of the simple geometry, one 
can perform the required integrations in closed form for each subclement and add 
the contributions. The details were described in Reference (7] for a four node ele- 
ment. It was shown in Reference [10] that the closed form expressions for the K {J 
arc quite simple for a four node element. 

Rectangular macro elements with rectangular subelements cannot accurately 
model wavy regions like that shown in Figure 4. For such microstructure one 
needs to use distorted subelements. In the more general case, such as when the 
interface between woven mats is not straight, the macro element will also be dis- 
torted. Figure 5 shows a distorted quadrilateral macro element with distorted 
subelements. The large numbers (1-4) are the macro element node numbers. The 
smaller numbers are the subelement node and element numbers. For simplicity 
the resin pockets are not modeled. 

To obtain a single field approximation, the subelement degrees of freedom 
(dof) must be expressed in terms of the macro element dof. There are several 
ways in which we can proceed. Two procedures will be discussed herein. Before 
proceeding it slum Id he (Jointed out that in general the single field character is 
only exactly satisfied at the sulteletucnt nodes. The first procedure is to consider 
the subclcincni mesh to be an ordinary finite clement mesh. Hie only difference 
is that after the subclcincni stillness matrix and equivalent nodal load vector are 
determined, they arc not immediately assembled, but are first transformed. This 
transformation can Ik expressed in matrix notation as 


3 



(4) 


K., = K'„. T nj 

r. = T„.r„ 


where 7._. is defined by q‘ = T,„q„ and 

q\ = nodal displacements for subelcmcnt 
q„ — nodal displacements for macro element 
= stiffness matrix for subelcmcnt 

K.j = subelemcnt contribution to stiffness matrix for macro element 


The transformation matrix T,„ is calculated using the macro element interpola- 
tion functions (which are defined in terms of local coordinates f and 77 ) evaluated 
at the subelement nodes. For example, for a four-node macro element and a 
three-node subelement the transformation is 


u\ 

v', 

u’l 

V‘l 


where t tJ = 


f» 

hi 

hi 


til 

hi 

til 


Njtft.nt) 0 

0 


til 

hi 

til 


ti 4 

h* 

1 34 



1 


v, 




v 4 


(5) 


Another possibility involves transforming the interpolation functions. This alter- 
native is much more efficient unless there are a very large number of integration 
points. This procedure will be illustrated by considering the interpolation for the 
displacement in the x-direction, u. A few more definitions are required before 
proceeding. 


u = macro element displacement in jr-direction 

u, = macro element nodal displacements in jr-direction 

u ’ = subelcmcnt displacement in jr -direction 

u\ = subelcmcnt nodal displacements in x -direction 

N. = interpolation functions for macro clement 

<v; = interpolation functions for subelemcnt 

Within ;i subelcmcnt the v -displacement is approximated as 
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«' = N\u\ 


(6) 


But the subclcmcnt nodal displacements arc slaves to the macro element nodal 
displacements, as described earlier. This can be expressed as 

«: = (7) 


where s'., 17, — coordinates of subclcmcnt node /. Combining Equations (6) and 
(7) gives 


!?.-)«> ( 8 ) 

or 

u' = N'.TijU; (9) 

where = Nj($,,T)i). Note that this transformation matrix T 0 is similar to that 
in Equation (4). The approximation for u can also be expressed in terms of 
modified interpolation functions, 


tt‘ = NjUj 


00 ) 


where Nj = N‘ T ;j . 

Since the range of i in Equation (10) is 1 — (number of nodes in the subele- 
ment) and the range of j is 1 — (number of nodes in the macro element), the 
“modified” interpolation functions can be different in number than the original 
functions. These modified interpolation functions are used when calculating the 
subelement stiffness matrices. Recall that the B matrix contains derivatives of the 
interpolation functions Nj. This presents no problem since the T tJ contains only 
constants. For example. 


dx dx " 


(ID 


These modified interpolation functions are used in evaluating the terms related 
to the displacement interpolation. The unmodified interpolation functions are 
used to determine the determinant of the Jacobian for use in mapping the differ- 
ential area d $dn from the subclement local coordinate system to a global coor- 
dinate system. Since the subclcmcnt displacements are slaved to the macro ele- 
ment displacements, there is considerable freedom in defining the subclcmcnts. 
For example, there is no need to prevent “dangling” nodes like that shown in 
Figure 5. In fact, one can even define the stillness matrix tor a macro clement to 
be a summation of some very unlikely looking subclcmcnts. This is shown sche- 
matically in Figure 6. 1 Ins is probably ol little practical utility for nvo- 
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dimensional models, but for three-dimensional models this represents a major 

The single field approximation gives very poor results for some configurations 
For example, if the lamina in Figure 3 have large differences in E r , it is very 
difficult to approximate the stiffness in the y-direction using a single field approx- 
imation. This is because the single field assumption results in continuity of 
strains, which causes a discontinuity of stresses which should be continuous at 
the lamina interfaces. A numerical example of this poor performance will be 
given in the “Results and Discussion” section. However, as will be illustrated 
later, there are realistic configurations with significant inhomogeneity for which 
a single field approximation performs well. Also, the macro elements described 
herein cannot be evaluated using the usual mesh refinement convergence meth- 
ods. As the mesh becomes more refined, the inhomogeneity within an element 
disappears and the macro element becomes an ordinary element. 

RESULTS AND DISCUSSION 

Results for two basic configurations will be presented. The first is a one- 
dimensional bimaterial rod and the second is a 2D idealization of a woven textile 
The material properties for the woven textile were assumed to be 

= I X) GPa Eii = 10 GPa E i3 = 10 GPa 

i' l2 = 0 35 v l} = 0.35 j/,, = 0.3 

G ti = 5 GPa G,i = 5 GPa G ls = 3.845 GPa 

These material properties are meant to represent those for a transversely iso- 
tropic tow. They do not correspond to any particular material system Two- 
dimensional material properties were obtained by imposing plane strain condi- 
tions The maieri.il properties were transformed to account for the inclination of 
the fiber bundle 

The bimaterial rod (shown schematically in Figure 7) was used to evaluate the 
accuracy of a single field approximation when two materials are loaded in series. 
The axial displacement was assumed to vary as T.% l a.x\ where n equals the 
order of the polynomial. Figure 7 shows the error in predicted stiffness verus the 
ratio EJE .. As expected, the error increases with the ratio EJE.. Perhaps sur- 
prising is the inability of an eighth order polynomial to adequately predict the re- 
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sjxmsc when /:»//:. is larger Ilian alKHit 2. Olwiously. ilie single field approx ima- 
lion is not very useful when l wo very dilTcrcnl materials arc loaded in series. 
However, most realistic configurations involving dissimilar materials have load 
paths which are a combination of series and parallel. The example of primary 
concern in this paper is a textile composite, which will lie discussed next. 

Two -dimensional idealizations of textile composites were anal y/.cd using single 
field macro elements. The low path was assumed to Ik sinusoidal. The thickness 
of the tow, /;/ 2, was kept constant along the path. Wavincss ratios bla (see sketch 
in Fiuurc 8) were varied from .083 to .333. It should Ik noted that a woven com- 
posite is inherently three-dimensional. There is no typical cross section. Con- 
comitantly, results from any two-dimensional textile model must be used with 
caution. Consequently, the results presented should only be interpreted as an 
evaluation of the effectiveness of the macro elements for handling microstruclurc. 
Figure 8 shows the variation of cxtcnsional stiffness with wavincss. Two symmet- 
rically stacked mats were considered. Only one mat was modeled. Symmetry 
conditions were imposed on the lower surface of the mat. Results were obtained 
using 60 eight-node traditional finite elements (reference solution) and 2 eight- 
node macro elements. The macro elements predict the stiffness variation quite 
well, except for very large waviness ratios. 

Figure 9 shows undeformed and deformed finite element meshes for a single 
textile mat using 8-node traditional and 12-node macro elements. This configura- 
tion is different from that in Figure 8, which had symmetry on the lower surface 
of the mat. The absence of symmetry constraints results in large bending defor- 
mation. The deformed meshes are also shown overlaid to compare the predicted 
shapes. The macro elements predict the deformed shape very well. 

Figures 8 and 9 showed the good performance of the macro element for pre- 
dicting global response. This does not imply that stresses or strains within the el- 
ement can be calculated accurately. In fact, the errors can be quite large. Figures 
10 and 11 show the variations of a t along the lower boundary of the axial tow for 
two symmetrically stacked mats. Results are included for both traditional and 
macroelement analyses. The sample points are labeled in the figures as points 
1-6. Figure 10 shows a x for a waviness ratio of .333. The actual o t variations, i.e. , 
that calculated using conventional finite elements) is not complicated, but the 
single-field approximation is quite inaccurate. A waviness ratio of .333 is fairly 
large. For a smaller wavincss ratio of .166 (Figure 11) the accuracy of the single- 

licld approximation is much l>ctter. However, the use ol single-field finite ele- 
ments to calculate local stresses and strains is not recommended. Much belter es- 
timates lor local stresses and strains can be obtained using a global/local strategy. 
Single-field macro elements can Ik very useful for the global analysis. A refined 
traditional finite element analysis can then be used for the local analysis. 
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CONCLUSIONS 

A new type of (mile element was developed for analysis of textile composites. 
1 lus new element (referred to herein as a macro element) accounts for the spatial 
variation of material properties within a single element. Tests of the macro ele- 
ments showed good performance for modeling the global deformation behavior 
ot textile composites. Because of die single field assumption, the stresses 
calculated inside the macro element arc not accurate. To obtain these stresses a 
global/local strategy should be used in which macro elements arc used for the 
global analysis and conventional finite elements arc used for the local analysis 
Although only two-dimensional elements were evaluated, the formulation is 
valid for three dimensions. However, there are challenges in 3D modclino which 
arc not so apparent or do not exist for 2D models. For example, in 3D one could 
imagine mats which are oriented at other than 0° or 90° relative to the macro ele- 
ment axes. Such an off-axis mat is much more difficult to model, particularly if 
it is combined with mats with other orientations. There is obviously still much 
work required to develop a general textile composite analysis. 
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Chapter: 2 


MACRO FINITE ELEMENT USING SUBDOMAIN 

INTEGRATION 


SUMMARY 

For some heterogeneous materials, it is not practical to model the microstructure directly using traditional 
finite elements. Furthermore, it is not always accurate to use homogenized properties. Macro elements 
have been developed which permit microstructure within a single clement. These macro elements 
performed well in initial tests. 


INTRODUCTION 

In traditional 2D and 3D finite-element analysis, the material properties arc assumed to be 
constant or at least to vary smoothly within a single element. This is valid for most engineering 
applications because the microstructural scale (e.g. grain size in a metal) is very small 
compared with the element size. However, some materials exhibit a very coarse 
‘microstructure*, such as laminated or textile composites. Figure 1 shows a schematic diagram 
of a cross-section of a textile composite. Owing to the complicated geometry, textile composite 
structures are very difficult to analyse. To use the traditional finite-element method to solve 
this kind of problem, finite elements have to be defined such that the material properties vary 
smoothly in each finite element. This results in a very large number of elements. For laminated 
composites, the geometry is simpler, but if there are many laminae (for example, 50 laminae 
for a 6*25 mm laminate), modelling of individual layers becomes impractical because of the 
large number of elements. 

Material homogenization is one way of treating inhomogeneous materials. In this procedure, 
the spatially variable material properties are replaced by some ‘effective* homogeneous 
properties. The effectiveness of material homogenization, however, depends on the problem 
to be analysed. Material homogenization theories 1 assume that the applied loading on the 
boundary of the representative volume element (RVE) is spatially homogeneous. This 
assumption is good as long as the characteristic scale of the microstructure is much smaller 
than that of the macrostructure. For example, volume averaging in laminated composites 
works well for in-plane loads. However, it gives large errors for bending loads unless there are 
many plies and the different plies are highly dispersed through the thickness. Higher-order 
theories such as classical laminate theory (CUT) 4 account for the geometric details of the 
microstructure for the laminated composite plates in terms of 1st and 2nd area moments of 
inertia. CLT works well for in-plane and bending problems in thin laminated plates. But CLT 
ignores out-of-plane strains, which is unacceptable for relatively thick plates. There are many 
other ways of homogenization, but none of them are problem-independent. When there are 
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material property discontinuities, it is most accurate to model each material group discretely. 
However, this approach requires large computer memory and CPU resources. 

Some work has been done in dealing with specific problems to overcome this difficulty. 
Steven 5 developed a quadratic triangular and a quadratic isoparametric element with an 
internal interface modelled by a straight line. In his work, trian ;ular subregions were used to 
perform numerical integration. He also suggested the possibility of using a second 
isoparametric mapping to simplify the integration, but he neither described any details nor 
implemented the method. The concept of splitting the integration limit was also discussed by 
Panda et a/. 4 In his finite-element formulation for laminated plates, the integration limit 
through the thickness was divided to define material properties of each individual layer. Foye 7 
studied material properties of fabric-reinforced composites using subcell analysis. The unit cell 
was divided into rectangular paraleltepiped subcells. Since the subcell boundaries do not match 
the material interfaces, averaged material properties were used in each subcell. 

The present paper describes a ‘macro’ finite element which can account for the details of 
microstructure within an element. A macroelement is defined to be an element consisting of 
several subdomains. Figure 2 shows a macro element that has four subdomains. Macro 
elements can have material discontinuities inside the element, but in each subdomain the 
material properties are smooth functions of the spatial co-ordinates. The macro element is 
identical to a traditional finite element when it has only one subdomain. When there are 
material discontinuities, the subdomainss are used to define the material boundaries and to 
facilitate the numerical integrations. 

It should be noted that since the present study was based on a displacement formulation, 
even with very high-order interpolation, significant stress errors are expected near the region 
where geometric or material discontinuities occur. Local stress distributions for regions of 
special interest can be achieved by global/local analysis. 4 " 10 The proposed macro elements are 
best suited for use in the global analysis. 

In the following Sections, the finite-element stiffness matrix formulation is explained in 
detail. Then examples for several configurations are discussed to illustrate the performance of 
the macro elements. 


CALCULATION OF MACRO ELEMENT STIFFNESS MATRIX 

This section describes the finite-element stiffness matrix formulation for a macro element for 
two-dimensional elasticity analysis. The extension to three dimensions is trivial and will be 
discussed briefly at the end of this Section. 

In a traditional displacement -based finite-element method, the dement stiffness matrix has 
the form 

(K)=[ (B) T (D) (B| dx dy (I) 

J o 
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where (D] and (B] are defined by 

lo| * IDHcJ 
Ul « (B](q| 

and (q| is (he nodal displacement vector. Supposing that there are material property 
inhomogeneities within the integration domain fi (i.e. the (D] matrix is a piece-wise-continuous 
function of spatial co-ordinates within a macro element), the macro element is divided into 
subdomains. Within each subdomain, the material properties vary smoothly (see Figure 2). 
Consider an element that has n subdomains Q, where 

S a, = o (2) 

»m | 

The element stiffness matrix becomes 

[K]=E( dxdy (3) 

A robust procedure is needed to evaluate the contribution of the fth subdomain 


IK], = l |B] T I O],- 1 B] dx dy (4) 

Jo, 

The procedure developed herein involves the use of three co-ordinate systems. Figure 3 shows 
the three co-ordinate systems. The use of three co-ordinate systems differs from conventional 
finite elements, which use a global and a local co-ordinate system. (In Figure 3, these are the 
(x,y) and ({, i|) systems.) The mapping of co-ordinate systems in conventional finite elements 
permits integration over a simple square region even when the actual finite element is quite 
distorted. If the material properties vary discontinuously within an clement, subdomains must 
be defined in which the material properties vary continuously. In general, these subdomains 
are distorted such as the one indicated by the shaded region ( ijkl ) in Figure 3(a). When this 
distorted subdomain is mapped into the ($, iy) co-ordinate system, it is still distorted (shaded 
region in Figure 3(b)) and concomitantly the integrations would not be simple. If this 
subdomain is mapped again into a third co-ordinate system ( r 9 s ), the integrations are again 
quite simple. The remaining task is to describe how to perform the integrations in the (r,s) 
co-ordinate system. 

There are two primary concerns. Tbe first is defining the differentia] element dx d y in terms 
of dr d 5. Figure 3 shows that 

dxdy — | J 1 d£ dij (5) 

where J is the Jacobian matrix defined by 


However, 


where 


3(x,y) 


d* dij= | J| drds 


( 6 ) 


d(r,s) 
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Therefore, the net result is 


dx dy = | J 1 1 j | dr d? (7) 

The second concern is defining the integrand in terms of the (r.s) co-ordinates. That is, 
fB] [DMB] involves derivatives of the interpolation functions. These interpolation functions 
arc defined in terms of £ and ij, not r and s. For example. 


’dNi' 


~ 9 N,~ 

9x 

= J-‘ 


BNi 


9Nj 

dx 


_ 


(8) 


Note that the calculation of the derivatives involves J but not J. This is because the .V,’s are 
defined in terms of £ and 7 alone. It is necessary to evaluate the integrand £B) T [D)«- [B] at 
particular values of r and s. As part of the mapping between the (£, ,) and (r.s) co-ordinate 
systems, £ and jj are approximated as 


* = N.(r.s)£, 

*1 = M(r,s)ij. 


When performing numerical integration in the (r.s) co-ordinate systems, £ and f are 
determined using equation (9). These values are then used to evaluate the integrand. 

With equations (4) and (7), the contribution of the rth subdomain stiffness matrix becomes 

IKJr* J' t J‘ [B] T (D]r(B] | J 1 1 J | dr ds (10) 

In Figure 2, both the macro dement and the subdomains are quadrilaterals. This is not 
necessary. The interpolation for the solution is defined in the ({, ij) co-ordinate system. The 
subdomain, which is mapped into the ( r 9 s ) system, is needed only to simplify the numerical 
integrations. The type of subdomain does not affect the solution except that it should 
adequately describe the geometry of the microstructure. 

The extension to three dimensions is simple. The three co-ordinate systems would be 
(x t y % z), ({, ij, f) and (r # $, (). The form of equations (1)— (10) is unchanged except to account 
for an additional co-ordinate direction. The contribution of the rth subdomain to the macro 
element stiffness matrix would be 


[K], = 



(B] T ( D], [B] | J | | J ( dr ds d/ 


(ID 


RESULTS AND DISCUSSIONS 

This section discusses the use and performance of the macro elements for two-dimensional 
elasticity. Three basic configurations were studied: (1) square and distorted 4-node elements, 

(2) (O/9O2/O) and (90/02/90) laminated beams with end moment and shear force loadings, and 

(3) a single and double plain weave textile composite under tension. Four-, 8- and 12-node 
macro elements were evaluated. The following material properties were used: 

£11 = 100 GPa £22 = 10 GPa £33= 10 GPa 

*12 = 0*35 *13 — O’ 35 *23 = 0*3 

C12 = 5 GPa G13 = 5 GPa G23 = 3 * 845 GPa 
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Two-dimensional material properties were obtained by imposing plane strain conditions. For 
textile composites, the material properties were transformed to account for the inclination of 
the fibre bundle. 

The first use of the macro element was to demonstrate that the mapping is correct. Square 
and distorted 4-node macro elements were subdivided into four subdomains (Figure 4). All 
four subdomains were assigned the same material properties. The macro element stifTness 
matrices should be the same as that for 4-node traditional elements. Tables I and U list the 
eigenvalues of the traditional and macro clement stifTness matrices for different orders of 
integration. For the square element shown in Figure 4(a), both traditional and macro elements 
produce exactly the same results when (2x2) Gaussian integration is used. Note that since 
there are four subdomains, the actual number of integration points for a macro element is four 
times the number of subdomain integration points. For the distorted element shown in 
Figure 4(b), results for the traditional finite element using (2 x 2) integration differ from the 
exact solutions. Table II shows that (3 x 3) integration for the traditional element and (2 x 2) 
integration for the macro element are nearly exact. As the integration points increase, both 
elements converge to the same results. 

Figure 5 shows the moment resultants for two laminated beams. Tip displacements were 
applied to produce a maximum strain of 0-14%. Homogenized material properties were 
obtained by the rule of mixtures. 4 The reference solutions are from traditional finite-element 
analysis with four 8-node traditional dements. Four-node and 8-node macro demenu were 
evaluated. A single macro element with four subdomains (one subdomain per lamina) was 
used. For the 4-node macro dement, selective ‘reduced* integrations 1 1 with 17 integration 
points were performed. ResulU show that for both (0/90j/0) and (90/0j/90) stacking sequences, 
one 4-node or 8-node macro dement predi. ts the bending stiffness very wdl. Of course, an 
8-node traditional dement with the volume-averaged homogenized material properties cannot 
distinguish differences in the stacking sequence. Hence, the errors are large for the volume- 
averaged homogenized material model, as expected. The percentage errors are shown in 
Table III. Note that accuracy depends on stacking sequence. 

Figure 6 shows the tip displacement comparisons for a short (3 x 1) cantilever beam for two 
stacking sequences. A unit shear force was applied at the right end of the beam. Four 
subdomains were used to account for the inhomogeneous material properties. Single 4-nodc, 
8-node and 12-node macro elements were used. The reference solutions were obtained with a 
refined mesh (64 dght-node elements). As expected, the traditional finite-element analysis using 
the refined mesh with volume-averaged homogenized material properties docs not predict the 
deformation behaviour. The 8-node macro dement predicted the displacements fairly well. The 
12-node macro element showed excellent performance. The 4-node macro element did not 
perform well. This was expected since the assumed displacement fields for the 4-node element 
are too simple for this problem. Table IV shows the percent errors for each case. For all three 
macro elements, the error was larger for the (90/02/90) laminates. 

The ‘effective’ extensional modulus E x against waviness of plain weave textile composites 
was calculated using traditional and macro elements. Figure 7 shows the configuration studied: 
two symmetrically stacked layers. Thick and thin lines in the upper mat indicate the macro 
elements and subdomains, respectively. Only the upper mat was modelled because of 
symmetry. The models have a length which is the same as the fibre bundle wavelength. For 
simplicity, the textile composites were assumed not to have any pure matrix regions. 
Displacements were applied to produce a 0*1 percent nominal strain <c> in the x-direction. 
The effective E, was defined to be 


where 


Em * 


<£> 

<C> 


<a> = 


Axial force 
Area 
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Th« waviaess was defined to be w = bja, where the centreline of the wavy fibre bundle is 
assumed to have a sinusoidal shape given by 



Two 4-, 8- and 12-node macro elements were used. Each macro element consisted of 15 
subdomains and models a half-wavelength. For the reference solution, a mesh with 60 
traditional eight-node elements was used. Figure 8 shows several of the traditional finite- 
element models. (The wavy fibre bundles are indicated by the shaded region.) Figure 9 shows 
the effective E, against waviness. The error increased with increased waviness. Both 8-node 
and 12-node elements performed fairly well. The 4-node element was not very accurate except 
for small waviness. 

Deformed meshes for macro element and traditional models are shown in Figure 10. These 
models are for a single plain weave mat (i.e. no symmetry). The Figure shows quite graphically 
the effect of the microstructure oa the predicted deformation of a single mat. The Figure also 
shows that the macro element predicts the deformed shape very well. It should be noted that 
only linear analysis was performed in the present study. The deformation shown in Figure 10 
is larger than would be expected from a non-linear analysis. 

CONCLUSION 

A displacement-based macro dement was developed to expedite elasticity analysis of 
heterogeneous materials. Two-dimensional macro elements with four, eight and 12 nodes were 
implemented and evaluated for several realistic configurations. Since the macro dements used 
a continuous strain field approximation, it is obvious that there is violation of equilibrium at 
the material interfaces and the stress distributions near the interfaces would not be very 
accurate. However, the macro dements performed well in terms of global response for the 
configurations considered. To obtain detailed local stress distributions, a global/ local strategy 
is needed. The proposed macro dements should be very useful for expediting the global 
analysis. 
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Tabic I. Eigenvalues against Gauss integration 
points for square dements (Figure 4(a)) 


Eigenvalues (x 10*) 


Integration 

1 x 1 

2x2 

4-node traditional: 

0 00000 

0-00000 


0*00000 

0-00000 


0*00000 

0-00000 


0*00000 

5*76923 


0*00000 

5*76923 


7*69230 

7-69230 


7*69230 

7-69230 


19*23076 

19-23076 

4-node macro: 

0*00000 

0*00000 


0-00000 

0-00000 


0*00000 

0-00000 


4*32692 

5*76923 


4*32692 

5*76923 


7*69230 

7*69230 


7*69230 

7*69230 


19*23076 

19*23076 


Table II. Eigenvalues against Gauss integration points for distorted elements (Figure 4(b)) 


Integration 


Eigenvalues (x 10*) 


I x I 

2x2 

3x3 

4x4 

5x5 

4-node traditional: 

0*00000 

0*00000 

0-00000 

0-00000 

0-00000 


0*00000 

0*00000 

0-00000 

0*00000 

0-00000 


0*00000 

0*00000 

0*00000 

0-00000 

0-00000 


0*00000 

4*55397 

4*61494 

4-61621 

4-61624 


0*00000 

6*19357 

6*25265 

6*25377 

6-25379 


7*18992 

8*56256 

8*62312 

8-62441 

8-62444 


7*93269 

9*38271 

9*46246 

9-46405 

9*46409 


20*57450 

20*87012 

20*87916 

10*87935 

20*87935 

4-node macro: 

0-00000 

0-00000 

0-00000 

0-00000 

0-00000 


0-00000 

0*00000 

0*00000 

0*00000 

0*00000 


0-00000 

0*00000 

0-00000 

0-00000 

0*00000 


7*93509 

4*61146 

4 61621 

4 61624 

4-61624 


8*59403 

6*24930 

6*25377 

6-25379 

6*25379 


3*63243 

8*61957 

8*62441 

8-62444 

8-62444 


5*06974 

9*45776 

9*46405 

9-46409 

9-46409 


20*76932 

20*87863 

20*87935 

20-87935 

20-87935 
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Tabic III. Percentage error for moment 
resultants 



(0/90j/0) 

(90/0,/ 90) 

4 -node macro 

2-451 

5*832 

8 -node macro 

0-001 

0-000 

12-node macro 

_ 


Homogenized 

37*61 

160*3 


Table IV. Percentage error 
for tip displacements 


(0/90^0) 

(90/0^90) 

23*78 

29*45 

9*873 

12-87 

2*857 

4*777 

33*11 

54*64 
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Figure I. Schematic diagram of 



woven composite 



X 

Figure 2. Typical 4-node macro element 
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4 8 12 16 0 4 8 

Moment, N-m Moment, N-m 

(t) (0/90/90/0) (b) (90/0/0/90) 

Figure 5. Moment resultants for two stacking sequences 
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Displacement, m Displacement, m 

(a) (090/90/0) <b) (90/0/0/90) 

Figure 6 Tip displacements of (3 x 1) cantilever beams of two stacking sequences 
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wWomain macro elements 


/ / \ 



Figure 7. Two symmetrically sucked plain weave mats 



Figure 8. Several traditional meshes: (a) A/e = 0 333. (b) A/e = 0 167. (c) 6/e = 0 1 II. (d) bio = 0 083 
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ooooo 2 four-node macro 

ooQoe 2 eight-node macro 

A 666 6 2 twelve-node macro 

$ 0 0 0 0 60 eight-node traditional 


0 
0.05 


0.10 


0.15 


0.20 


0.25 


0.30 


0.35 


Waviness, b/a 

Figure 9. Exiensional modulus against waviness 


(a) Traditional mesh 


<b) Macro element mesh 



(c) Overlaid deformed meshes 

Figure 10. Deformed meshes for a single plain weave mat (120 eight -node traditional elements and 4 twelve-node 
macro elements are oscd; waviness bja =0*167: nominal strain =0*05) 
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Chapter: 3 


ENHANCED DIRECT STIFFNESS METHOD FOR 
FINITE ELEMENT ANALYSIS OF TEXTILE COMPOSITES 


ABSTRACT 

Traditional homogenization techniques are not useful when the microstructural scale of 
a material is of the same order of magnitude as the structural scale of a component. Such is the 
case for many textile composites. Since discrete modeling of the microstructure throughout a 
component is prohibitively expensive, continuum finite elements are needed which account for 
microstructure within a single element. This paper describes a simple substructuring technique 
for formulating these special elements. 


INTRODUCTION 

By changing stacking sequence, fiber orientation, and materials, traditional composite 
laminates can be tailored for specific applications. With the introduction of advanced textile 
composites, there are even greater opportunities to tailor composite properties. Not only are 
there many textile forms (eg. weaves, braids, knits, etc.), but there are many unique varieties 
of each form. 


Accurate predictive analyses are essential for designing high performance composites. 
In contrast to traditional tape laminates, verified analyses are not in place tor textiles. Figure 1 
illustrates the complexity of the task of developing an accurate textile analysis. The figure shows 
schematics of a traditional laminate and a woven material. For the traditional laminate one can 
define a unit cell of dimensions approximately .007 mm. This unit cell can be analyzed to 
determine effective engineering properties for the much larger individual lamina. Then each 
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lamina can be treated as a homogeneous orthotropic layer. For a woven composite the unit cell 
can be larger than l mm. For the woven composite one can use homogenized engineering moduli 
to describe the tow properties (a tow contains on the order of 6000 filaments), but there is a 
much larger microstructural scale related to the interlocking of the tows. 

Description of the material properties for a weave requires different strategies at different 
microstructural scales. Figure 2 illustrates different microstructural scales. Actually all of the 
schematics were generated from the same basic unit cell labeled "coarse microstructure." The 
term "coarse" refers to the very distinct phases at this level of observation. In contrast, if a very 
large number of unit cells are considered, the material appears almost homogeneous (schematic 
labeled fine microstructure). At the extremes of microstructural scale the choices for material 
modeling are obvious. For coarse micro structure the individual tows and matrix pockets are 
modeled discretely both in terms of geometry and the abrupt changes in properties at the 
interfaces. For fine microstructure effective homogenized engineering properties can be used. 
Traditional finite element methods are appropriate at these two scales. Between these two 
extremes (labeled "transitional microstructure”) traditional finite elements are not appropriate. 

In this range there are too many microstructural features to model them all discretely, but there 
are too few to use homogenized material properties. In the transitional range of microstructure, 
special finite elements are needed which permit material variation within an element. Of course, 
this is routine for layered plate and shell elements. 

Recently, continuum elements have been developed for accounting for textile type 
microstructure within a single element [1,2]. The elements described in these references are 
based on a single assumed displacement field throughout the entire element. A more general 
element formulation is presented herein that includes the single field approximation as a 
degenerate case. This more general formulation is an example of reduced substructuring [3]. In 
brief, the implementation begins with the development of an ordinary finite mesh for the basic 
textile unit cell. Then interior degrees of freedom are statically condensed out. Next the number 
and location of desired boundary degrees of freedom are selected. Finally, the original boundary 
degrees of freedom are expressed in terms of the desired boundary degrees of freedom. One 
objective of this paper is to describe a very simple technique for calculating the stiffness matrix 
for a reduced substructure. The other objective is to show a few results which illustrate the 
effectiveness of this type of element. 

In the discussion that follows, the term "macro element" will be used to indicate an 
element which allows for internal microstructure. Accordingly, the elements described in [1,2] 
are single-field macro elements. Similarly, the reduced substructure elements will be referred 
to as multi-field macro elements, since the displacement field inside the macro element is defined 
piecewise. 


THEORY FOR REDUCED SUBSTRUCTURING 

In multi-field elements the internal dof are eliminated using the equivalent ot static 
condensation. Also, boundary degrees of freedom (dof) which are not to be part ot the macro 
element dof arc expressed in terms of tiie substructure dof using multipoint constraints. 
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Theoretically, this is all very simple. Con ider the 4-node macro element in Figure 3. 
Assume the governing equations are partitioned as follows 


k aa 



1a 


F a 

*Ib 

K bb 


1b 


f b 


and q A is the list of unknowns to be condensed out (see Figure 3). 


Before imposing the multipoint constraints on the excess boundary dof, the reduced 
stiffness matrix and load vector can be expressed as 


K,b - K bb ' *Ib Zl K ab 


F b = F b~ K ab K'a\ f a 


( 2 ) 


This procedure often is not practical as stated because of the matrix inversion which 
eliminates sparsity in K AA and the large matrix multiplications. The elimination of internal dof 
can also be accomplished using Gaussian elimination if the dof to be eliminated are grouped 
together at either the beginning or the end of the list of unknowns. This procedure is well 
known, so it will not be discussed herein. See [4] for details. After eliminating the interior dof, 
multipoint constraints can be applied to the remaining dof to eliminate unwanted boundary dof. 
This can be expressed in matrix form (assuming the four node macro element in Figure 3) as 


a B = T a where a 

“ macro “ macro 


‘U 


r ll 


w 


14 


(3) 


The transformation matrix T expresses how the excess boundary dof are slaved to the 
macro element dof. It should be noted that if the internal dof are also slaved to the macro 
element dof (rather than statically condensed), a single-field approximation is obtained. Of 
course, a formulation like that in [1] is much more efficient for single-field elements. However, 
the current formulation permits great flexibility for evaluating various approximations. 

It is not always efficient to order the dof such that Gaussian elimination can be used to 
obtain the reduced stiffness matrix and load vector, since such ordering might result in large 
bandwidths. An alternative is to use the formal definition of the stiffness coefficient K-. 

Kjj = force at dof i due to unit displacement at dof j 
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Using this definition we would simply solve a series of problems in which one dof is set equal 
to 1 and the rest of the boundary dof would be constrained to zero. The restraint forces at all 
the boundary dof constitute one column of the reduced stiffness matrix. 

This process is repeated for each boundary dof to obtain the entire reduced stiffness 
matrix. The reduced load vector is obtained by solving one additional problem in which all 
boundary dof are constrained to be zero and the internal loads are applied. The negative of the 
boundary restraint forces constitute the reduced load vector contribution for the internal loads. 
Once the reduced set of equation is obtained, the multipoint constraints can be imposed to 
eliminate unwanted boundary dof. This alternative is not new. It can be considered a numerical 
application of the direct stiffness method for calculating stiffness matrices. It also may not be 
very efficient when there are a large number of boundary dof to be eliminated. Consider a case 
in which there are 32 boundary dof, but only 8 are to be retained in the macro element. The 
procedure described above requires the solution of 32 unit displacement cases. A new procedure 
is discussed next which would only require 8 unit displacement solutions. 


ENHANCED DIRECT STIFFNESS METHOD 

The enhanced direct stiffness method is derived starting with a consideration of the work 
performed by the boundary nodal forces during deformation. To simplify the discussion only 
linear configurations will be considered herein. Figure 3 show a schematic of a typical mesh for 
a macro element. There are four interior nodes (nodes 1, 2, 3, 4), four boundary nodes to be 
retained_(nodes 11, 12, 13, 14; dof = qj), and six boundary nodes (nodes 5, 6, 7, 8, 9, 10; 
dof = q@) which are slaved to the qj through multi-point constraints. The nodal forces 
corresponding to q; and q^ are defined to be F; and F^, respectively. For the particular mesh 
in Figure 3, the range of i and /3 axe 1-8 and 1-12, respectively. Assuming linear elasticity, the 
work performed by the boundary nodal forces is 

* fa) • 

( 4 ) 

i * 1 , number of retained dof 
P = l, number of slaved dof 


The qp are slaved to the q v which can be expressed as 


(5) 


where T^-, = Nj(f^, n^) is calculated using interpolation functions for the boundary. Combining 
equations 4 and 5 yields 


W = 



( 6 ) 


It is well known that the stiffness matrix can be expressed as [5] 
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K 


dq m dq 


(7) 


Bui U = W for linear configurations so that 

K — 

" ' dq m dq a 


(S) 


Combining eqns. 6 and 8 yields 


dF 




dF 

m 

d 4n 


El 


'Pc 


dF. 

— ~ T * 

8 <?„ P " 


(9) 


8F 8F 

From the Maxwell-Betti reciprocity theorem we know that — - = — - . The third and fourth 

3 4,„ d <J„ 

terms in equation 9 are also equal for the same reason, but it is not obvious in the present form. 
To make the equivalence more obvious, first, equate the work of the forces F ^ with that of die 
equivalent forces f { in terms of the retained dof 


~ ■ q ' 


( 10 ) 


Combine equations 5 and 10 to obtain 




(ID 


Equation 1 1 shows that the equivalent nodal forces f - t are 


/. = 


( 12 ) 


Hence, the third and fourth terms in equation 9 become 
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Again, using the Maxwell-Betti theorem, these two terms are equal. Therefore, equation 9 can 
be expressed as 


K = 


dF. 


*f. 


d 4 m 


( 14 ) 


The implementation of equation 14 using the direct stiffness method is as follows: 

1 . Impose a unit displacement q j . 

2. Also impose displacements qp. since = T^jqj. 

3. Analyze model. 

4. Calculate restraint forces F n and F^. 

5. Calculate f n = F^T^ n 

6. The sum of F n and f n = column 1 of the reduced stiffness matrix. 

7. Repeat steps 1-6 for eacli q-,. 


CONFIGURATIONS 

Plain weave composites with different waviness were analyzed. Figure 4 shows a 
conventional 3D finite element model of a plain weave. It has 381 nodes and 64 quadratic 
elements. The tow path was assumed to the sinusoidal. The waviness ratio is defined to be b/a, 
where b = the mat thickness and a = the wavelength for the tows. The waviness ratio was 
varied from .033 to .33. 

This mesh was used to obtain reference solutions. It was also used to generate 20-node 
single-field and multi-field macro elements. Hence, there were three models: the conventional 
model shown in Figure 4, a one element mesh using a 20-node single field element, and a one 
element mesh using a 20-node multi-field element. The single-field results were obtained using 
the formulation in [1]. 

Two sets of boundary conditions were used: one for a narrow two mat composite and the 
other for an infinitely repeating unit cell. The boundary conditions for the narrow two mat case 
correspond to a specimen which is infinitely long in thex-direction, width "a" in the v-direction, 
and thickness 2b in the z-direction. The boundary conditions were as follows. 
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Narrow two mat composite: 


«(0,y.z) =0 

v(x,0.z)=0 
w(x,y,0) = 0 



= specified constant value 


(15) 


Infinitely repeating unit cell: 

Constraints listed in equation (15) and 

a ) 

, — ,z =constant 

2 J 

b) = constant 



The material properties for the tows and resin pockets were assumed to be 

Tows: 


E n = 206.9GPa 

"12 = -25 
G 12 = 2.386GPa 


E 22 = 5.171GPa 
"13 « -25 
G 13 = 2.386GPa 


Resin: 


E 33 = 5. 171GPa 

"23 = -25 
G 23 = 2.386GPa 


E = 3.45GPa v = .35 


RESULTS AND DISCUSSION 


There are two aspects to the evaluation of the procedures outlined in this paper. First, 
the methodology for calculating the multi-field stiffness matrix was checked. This was 
accomplished by comparing the stiffness matrix with that obtained using standard Gaussian 
elimination followed by application of multipoint constraints. As expected, the results agreed. 
The second task is to evaluate the performance of the multifield elements for analysis of textile 
composites. This second task is only partially complete. A few results are discussed in this 
section which suggest that this type of element can be very useful. 

Axial loading along the x -direction of a narrow strip of plain weave composite was 
modeled, as described in the Configuration section. Because of the complex spatial variation of 
materials properties, there is significant distortion, even under simple extension. Figure 5 shows 
the distortion of the macro element mesh and the conventional mesh. The macro element predicts 
the distortion quite well. ( It should be noted that the elements in Figure 5 are drawn with 
straight lines joining the nodes. This is a limitation of the plotting software, not a characteristic 
of the solution.) 
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Figures 6 and 7 show the variation of several effective engineering properties with 
waviness ratio for infinitely repeating unit cells. Results are shown for conventional, single-field, 
and multi-field elements. Both types of macro elements predict the trends quite well. As 
expected, the performance of the multi-field elements is considerably better than that for the 
single-field elements. The accuracy of the multi-field elements is quite good except for very 
large waviness ratios. At small waviness ratios the single-field macro elements predict the in- 
planc behavior very well, but not the out-of-plane ( ie . i^). The single-field approximation 

imposes strain continuity throughout the element, which is not correct for heterogeneous regions. 
The error associated this approximation is more significant for out-of-plane properties than for 
in-plane properties. 


CONCLUSIONS 

A simple formulation for multi-field continuum finite elements with microstructure was 
developed. Initial tests showed very good performance in modeling the global response of a plain 
weave composite subjected to axial extension. Much more work is needed to fully evaluate the 
performance of these elements. Future work is needed ( and is planned) to evaluate the accuracy 
of these elements for much more complex loadings. Also, planned is an evaluation of the 
accuracy of the calculated stress fields within the elements. 
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Traditional laminate 



U 

.007 mm 


Textile mat 




Repeat length > 1. mm 


Fig. 1 Comparison of microstructural scales for traditional and textile composites. 
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Coarse Microstructure 



Fine Microstructure 

















Waviness Ratio = b/a 


Figure 4 Original finite element mesh for textile composites. 
(381 nodes, 64 elements) 
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Conventional 
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Figure 5 Deformation of multi-field macro element and conventional finite element models 


Modulus, GPa 




Poisson Ratio 



Waviness Ratio, b/a 


Figure 7 Poisson ratio versus waviness ratio for infinitely repeating plain 
weave textile composites. 
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Chapter. 4 

Boundary Effects in Woven Composites 


Abstract 

Two dimensional finite elements were used to study boundary effects in plain weave 
composite specimens subjected to extension, shear, and flexure loads. Effective extension, shear, 
and flexural moduli were found to be quite sensitive to specimen size. For extension and flexure 
loads stress distributions were affected by a free surface, but the free surface boundary effect 
did not appear to propagate very far into the interior. For shear load the boundary effect 
appeared to propagate much further into the interior. 

Key Words: textiles 

woven composites 
finite elements 
stress analysis 
boundary effects 


Introduction 

Fiber tows, each consisting of thousands of individual filaments, can be woven, braided, 
knitted, etc. to create complex fiber preforms. These preforms are then impregnated with a resin 
and cured to make textile composites. The interlacing of the fiber bundles provides many 
obstacles to damage growth. Accordingly, there is the potential for greatly improved resistance 
to impact damage growth. Unfortunately, there are also negative effects due to the fiber tow 
interlacing. The fiber tow curvature reduces the effective in-plane moduli. The curvature also 
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induces many local stress concentrations which can result in early diffuse damage initiation, 
particularly in the matrix. 'Hie fabrication process is not benign. For example, weaving involves 
much mechanical handling of unprotected fibers (i.c. fibers which arc not embedded in matrix). 
Stitching of textile preforms to increase delamination resistance has the side effects of breaking 
fibers and inducing local fiber curvature. Optimal design requires the capability to predict both 
the positive and negative effects of potential textile fiber architectures. Unfortunately, the 
complex fiber architecture is difficult to analyze. Accurate analysis requires accurate geometric 
representation and constituent properties, such as fiber and matrix properties and fiber volume 
percentage. For textile composites there is particular difficulty in determining the actual fiber 
tow geometry and developing a three-dimensional model which can be analyzed. There have 
only been a few attempts at detailed three dimensional analysis (eg. Refs. [1-3]). Even the 
accuracy of these models for local stress calculation is an open question because of the 
uncertainties in the input data (i.e. the approximation of tow geometry and other properties). 
Most of the analyses to date have been similar to laminate theory in level of approximation or 
detailed two dimensional (2D) or quasi-three-dimensional (Q3D) numerical analyses of a 
"representative" cross-section (eg. Refs. (4-7)). As the schematic in Figure 1 shows, there is no 
such “representative" cross-section, even for a plain weave composite. While such 2D or Q3D 
analyses are likely insufficient for accurate prediction of local stress states, they are useful for 
obtaining insight about the effects of fiber tow waviness on effective moduli and strengths. In 
fact, the results in this paper, which are based on 2D analyses, fall into this category. 

The analysis of textile composites is in its infancy as compared to laminated composites. 
There are many aspects of the behavior of these materials which have not even been examined, 
much less accurately described. The objective of this paper is to begin to address one question 
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about the behavior of plain weave composites: “How docs the presence of a boundary affect the 
stiffness and stress distribution in a representative unit cell?" rhe boundary surfaces referred to 
here arc those present due to finite thickness. Three nominally simple boundary conditions were 
considered herein: in-plane extension, transverse shear, and flexure. Configurations of different 
thicknesses were analyzed using 2D finite elements. The analyses were performed using 
conventional elements and multi-field macro elements (reference 8). Macro elements are defined 
to be elements which contain internal microstructure. The multi-field elements arc a form of 
reduced substructuring. The macro elements permitted analysis of quite large models without 
requiring huge amounts of computer memory and epu time. Of course, a few macro elements 
are not as accurate as using a huge collection of conventional elements. Accordingly, one 
additional objective of the paper is to evaluate the performance of macro elements for simple 
configurations. 

The following sections will begin with a discussion of the configurations studied. Then 
the results will be discussed. First effective extensional, shear, and flexural moduli will be 
discussed. Then the effects of boundaries on stress distributions will be discussed. 

Configurations 

The various configurations studied are all synthesized from a single basic unit cell. This 
unit cell will be discussed first. Then boundary conditions for infinite and finite configurations 
will be discussed. 

Unit Cclj 

The basic unit cell is shown in Figure 2 . The cell consists of tows running in the x- and 
z- directions. In reality there would also be pure matrix pockets, but these were filled with z- 
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direction lows in the model used. (Of course, in reality there is no typical cross section either, 
as discussed earlier.) The two dimensional approximation implies dial the x- direction low is a 
wavy "plate" and the z-direction tows are straight fiber bundles. Obviously these arc serious 
approximations, so the results presented are intended to be qualitative only. The centerline of 

the x-direction tows follows a wavy path described by the function For the results 

presented herein a = 1.5/3. The thickness of the tow as measured along a line normal to the tow 


centerline was held constant. It should be noted that the uni: cell selected assumes a symmetric 

stacking of the woven mats. There are an infinite number or other possibilities. 

Two sets of two material properties 

were used. They are 


Set I 



E n = 100 GPa 

= 10 GPa 

= 10 GPa 

i/j2 = 0.35 

v l3 = 0.35 

"23 = °- 3 

G 12 = 5 GPa 

G 13 = 5 GPa 

G n = 3.845 GPa 

Set II 



E u = 165.8 GPa 

E^ = 11.51 GPa 

E 3 3 = 11-51 GPa 

i/ 12 = 0.273 

i- 13 = 0.273 

1-23 = 0.33 

G l2 = 15.4 GPa 

G 13 = 15.4 GPa 

G.J 3 = 4.17 GPa 

These properties were transformed to account for the waviness 

of the x-direction tow. 

Plane strain conditions were imposed 

to obtain two dimensional properties. Two sets of 


properties were used. This is admittedly not optimal. The homogenization analyses were 
performed using Set l. The stress analysis results were obtained using Set ll. 
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Periodic Boundary Conditions (or Infinite Configurations 


Figure 2 shows a typical unit cell for symmetrically stacked mats before deformation. 
If this cell is imbedded within an infinite array of identical cells and displacements or tractions 
are imposed "at infinity", then every unit cell will deform identically. The periodicity of the 
displacement field can be imposed on a single unit cell, thus permitting the solution for the 
infinite domain. The solution for an infinite domain will be useful for comparison with finite 
configurations subjected to nominally uniform extension or shear. Using the coordinate system 


in Figure 2a, the periodic conditions can be expressed as 


u(or,y) = u(-a,y) + u 2 - Uj 

(1) 

v(a,y) = v(-a,y) + v 2 - Vj 

(2) 

u(x,/3) = u(x, -0) + u 4 - Uj 

(3) 

v(x,0) = v(x, -0) + v 4 - v t 

(4) 


There are no specified non-zero forces (The net forces are zero at any point inside the infinite 
media.). The "load" consists of the values chosen for (u 2 - u t ), (v 2 - v t ), etc. These values 
depend on the nominal strain state desired. (Specific values for the different states will be 
discussed later in this section. Equations 1-4 impose certain constraints which are not so 
obvious, but are worth mentioning, since they are exploited in the finite element analysis. These 
constraints are 


u 3 - u 4 = u 2 - u t 


v 3 - v 2 = v 4 - V, 


U 3 - u 2 = u 4 - U t 


v 3 - V 4 = v 2 - V, 


(5) 

( 6 ) 
(7) 
(S) 
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These constraints can be obtained from equations 1-4 by substituting in specific vertex values 
of x and y. Tor example, substitute x=cr into equation 3. 

u(a,0) = u(ar, -0) + u 4 - u, 

But u(a,0) = u 3 and u(a, -0) = u 2 . Hence, equation 3 slates that u 3 - u 2 = u 4 - u,. Equations 
5-8 indicate that if the nodal displacements at the four corners of the unit cell arc used to 
calculate the displacement gtadieuts, we find that (!(.(£(■(§( • - (£[ are COTSCUrt ' The 

subscript ”0” is used to indicate that these are nominal displacement gradients. On a pomtwise 
basis these are certainly not constant for the obviously inhomogeneous unit cells. Equations 1-4 
can now be expressed as 

u(a,y) = u (-a,y) + 3a 

v («.y) = v (-«.y) + 

u(x,0) = u (x, -0) + 20^^j 

v(x,0) = v(x, -0) + 20 ( 12 ) 

Because of symmetries only part of the unit cell must be modeled. Herein the quarter unit 
cell shown in Figure 2(b) was modeled. If all the symmetries had been exploited, only 
one-eighth of the unit cell would have to be modeled. For convenience the coordinate system 
is shifted to trie center in Figure 2(b). 
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For extension loading the boundary conditions are quite simple. The constraints imposed 


for nominal c x loading arc 




= specified constant value 


= constant, but unknown 


(13) 


Nominal a y loading (which was not considered herein) would be very similar. For 


nominal a_ v load the boundary conditions are 


“(*’ 2 ) ' -“(*'2) * 


= specified constant value 


v | — — , y | = -v | | = specified constant value 






(14) 


v (- x -ir) ■ -”R) 


The boundary conditions in equations 14 state that the displacements normal to an edge 
are anti-symmetric (and unknown except at the vertices). The tangential displacements are 
constant along an edge and are specified. 

Boundary Conditions for Finite Confintraiions 

Extension, shear, and flexure loading were considered for a wide range of specimen 
thickness (in the y -direction). Hence, the various meshes had different numbers of unit cells. 
For extension loads the boundary conditions were like those in equation 13 if one considers or 
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and (i to be the dimensions of the entire mesh, rather than just a quarter unit cell except that the 
top surface was traction free. Hence, the normal displacement “v* was not constrained to be 
constant along the top. For shear load all boundary displacements were constrained to follow the 
deformation u = cy and v = ex. Consequently, the boundaries remained straight after 

deformation for shear loading. 

For flexure loads the top and bottom surfaces of the model were traction free. A linear 
variation of normal displacements were imposed on left and right ends of the model. 

Results and Discussion 

There are two types of results which will be discussed. The first will illustrate the effect 
of specimen thickness on effective moduli. The second will illustrate the effect of unit cell 
location on stress distributions. 

E ffective Moduli 

For nominally simple deformation states, the effective engineering properties are expected 
to converge to constant values as the specimen thickness increases. Figure 3 shows the variation 
of the normalized effective E*. Figure 3a shows the variation of the average E* with the number 
of unit cells. The E x is normalized by the E x for an infinite array of unit cells modeled using 
conventional finite elements. The three curves were obtained using conventional finite elements 
and 8-node and 12-node multi-field macro elements. The 8-node macro element must be 
inherently a little too stiff, since it converges to a value approximately one percent too large. 
The 12-node macro element agrees very well with the conventional finite element results. For 
8 unit cells through the thickness the effective E x is within about one percent of convergence. 
Tins indicates that a specimen would need to be 8 unit cells thick to give an effective E x within 


49 



one |>erccm of a very ihick specimen. Figure 3b shows the variation of the effective E x with 
position for a configuration which has eight unit cells through the thickness. The effect. vc E x 
for each quarter unit cell was calculated based on the strain energy in the region. This is not a 
rigorous definition, but it docs offer some insight. The figure shows that the boundary quarter 
unit cell is about 18 percent softer than an interior quarter unit cell. The next quarter unit cell 
is about 5 percent too stiff. The third quarter unit cell has almost exactly the same stiffness as 
cells which are much further from the boundary. There is an obvious boundary effect, but it dies 
out very quickly. 

Figure 4 shows the effect of model size on normalized effective shear modulus G xv . 

In contrast to E x , the shear modulus converges from the stiff side. This difference is a 
consequence of the boundary conditions imposed. For E x there were free surfaces. The traction 
free condition permitted warping deformation to occur more easily near the free surface than in 
the interior, so the boundary caused softening. In contrast, all of the finite size shear specimens 
had specified x- and y- displacements over the entire boundary. This fully constrained boundary 
deformation resulted in larger effective G xy for smaller specimens. Figure 4 also shows that 
8-node macro elements perform poorly in shear. The 12-node macro elements perform quite 
well. It is interesting to note the distribution of the strain energy in a finite size shear model. 
The bar chart in Fig. 5 shows the strain energy in each quarter unit cell for a 3x3 array of unit 
cells. The effect of the boundary on the strain energy distribution is obviously quite complex. 

Figure 6 shows the variation of normalized flexural modulus with model size. The 
flexural modulus is defined to be (flexural stiffness)/!, where 1 = the second moment of the 
area. The flexural modulus in Figure 6 is normalized by the value for a configuration which is 
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ten cells thick. The flexural modulus converges more slowly than the cxtcnsional modulus, 'llic 
12-nodc macro element performs very well. The 8-nodc macro element is a little too stiff. 

Stress Distributions 

Figures 7-9 illustrate the effect of a free surface on stress distributions. Distributions are 
shown for extension, shear, and flexure. The stresses shown are evaluated with respect to the 
xy (global) coordinate system. 

Figure 7 shows the stress distributions for extension loading for three unit cells from two 
different configurations. One configuration had two unit cells through the thickness. The other 
had six unit cells through the thickness. The locations of the unit cells considered are indicated 
by shading in the figures. The waviness of the x-direction tow and the inhomogeneity causes a 
complicated variation of all three stresses. The <r x variation in the longitudinal tow is dominated 
by flexure induced by tow straightening, as shown by the locations of. maximum and minimum 
o x . The <r y is largest where the tows contact. The <r xy is largest where the tow rotation is largest. 

There are both striking similarities and differences in the stress distributions for the three 
unit cells. Figure 7 shows that the interior and exterior unit cells have very different stress 
distributions. There is .obviously a significant free surface effect. The exterior unit cells in 
Figure 7 have very similar distributions for. all three stress components. This suggests that for 
extension load the response of the exterior unit cells is not very sensitive to the total specimen 
thickness. 

The interior unit cell exhibits almost the same symmetries that one would expect from 
a cell embedded inside an infinite array. Also, the interior half of the exterior unit cells has 
stress distributions which are very close to those for the lower half of the interior unit cell. 
Apparently the free surface effect does not propagate very far into the interior. 
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Figure 8 shows ihc sircss distributions for shear loading. Single unit cell and 3x3 unit 
cell configurations were studied. Only the a y and o xy distributions arc shown, since o x was quite 
small. In this ease there are no free surfaces. (Displacements were specified along the entire 
boundary.) As was the case for extension, the interior and exterior response is different. The 
interior unit cell is located in the middle of the finite element model. Hence, the symmetries 
exhibited by the interior cell do not indicate the attenuation of boundary effects. In contrast to 
extension load, Figure 8 shows that for shear load the response of the boundary unit cells is very 
sensitive to total specimen size. Further studies are needed to determine the boundary layer 
thickness for shear loads. 

Figure 9 shows stress distributions for flexure loads. Only exterior unit cells are 
compared. The single unit cell model was subjected to a combination of extension and flexure 
so that the loading would be comparable to the exterior unit cell of the thicker model. The 
thicker model was subjected to pure flexure. Both models have free surfaces at both the top and 
bottom. The maximum c x does not occur at the free surface. This is because local flexure of the 
wavy fiber tow as it tries to straighten attenuates the o x . The top halves of the two unit cells in 
Figure 9 have very similar a x , a y , and a xy distributions. The lower halves exhibit much more 
differences. This is not surprising since the lower surface of the single cell is traction free but 
the lower surface of the cell from the thicker model is not. These results further indicate that 
there is a free surface effect (in this case, from the lower surface of the single unit cell model), 
but that the boundary layer is quite small. Finally, it should be noted that the stresses were lower 
for the flexure case than for the extension case even though the maximum nominal axial strain 
was .001 for both. 


52 



Conclusions 


Boundary effects were studied for woven composites subjected to in-plane extension, 
shear, and flexure. Effective moduli and stress distributions were calculated for configurations 
ranging from very thin to very thick. Only two dimensional models were studied. Since woven 
textiles are really three dimensional, these two dimensional results should only be interpreted 
qualitatively. Boundary effects were significant both in terms of stiffness and stresses. 

A specimen thickness of 6-8 unit cells was required to obtain moduli within about 2% of that 
for very thick specimens. For extension and flexure loading the stress distribution in exterior 
unit cells were quite insensitive to total specimen thickness. There appeared to be a characteristic 
response of boundary cells. Also, the boundary effect did not propagate very far into die 
interior. The response for shear load was more complex than for extension and flexure. Further 
work is needed to characterize boundary effects for shear loads. 
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a 

(b) Quarter unit cell 


Figure 2 Basic two-dimensional unit cell models. 
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Normalized extensional modulus E 



(a) Average normalized E x vs number of unit cells through thickness. 

Fig. 3 Normalized extensional modulus E x - Eight-node traditional elements 
were used for the infinitely repeating unit cell case. 
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Normalized extensional modulus 



(b) Normalized extensional modulus vs. position in an 8-unit 
cell configuration. (The sketch only shows four unit cells, 
since the configuration is symmetric.) 


Figure 3, completed. 
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Normalized shear modulus 



Number of unit cells through the thickness of the configuration 


Fig. 4 Normalized shear modulus vs. number of unit cells through the thickness 
of the configuration.(The number of unit cells is the same in both the 
x- and y- directions). 
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Figure 5 Normalized strain energy distribution in 3x3 unit cell model subjected to shear load. 

Strain energy in each quarter unit cell is normalized by that for an infinitely repeating 
unit cell array subjected to shear. 
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Normalized flexural modulus 



Number of unit cells through the thickness in the configuration 

Fig. 6 Normalized flexural modulus vs. number of unit cells through the 
thickness of the configuration. Results were normalized with the 
flexural modulus for a ten unit cell model. 
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(iii) Interior unit cell of model with six 
unit cells through thickness. 

(a) Axial Stress 


Figure 7 Stress contours for a two dimensional model of a plain weave composite 
under extension ( nominal axial strain = .001). 
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(i) Top unit cell of model with two 
unit cells through thickness. 



(ii) Exterior unit cell of model with six 
unit cells through thickness. 



(i ii) Interior unit cell of model with six 
unit cells through thickness. 


(b) Transverse Stress 

I'icurc 7, Continued. 
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(iii) Interior unit cell of model with six 
unit cells through thickness. 

(c) Shear Stress 


Figure 7. Concluded. 
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(iii) Interior Unit Cell of a (3x3) unit cell model 

(a) Transverse Stress 


Figure 8 Stress contours for a two dimensional model of a plain weave 
composite under shear, (nominal shear strain = .001) 
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(i) Single Unit Cell 


Stress, Pa 



4.000e+007 

3.666e+007 

3.332e+007 

2.998e+007 

2.664e+007 

2.330e+007 

l.995e+007 

1.661e+007 

l.327e+007 

9.932e+006 

6.591e+006 

3.250e-r006 



(iii) Interior Unit Cell of a (3x3) unit cell model 

(b) Shear Stress 


Figure S, Concluded. 
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(i) Single Unit Cell 




Stress, Pa 


L200c+008 

I.085e+008 

9.691e+007 

8.536e+007 

7.382e+007 

6.227e+007 

5.073e+007 

3.918e+007 

2.764e+007 

l.609e+007 

4.545e+006 

-7.000e+006 




(ii) Exterior Unit Cell of a (3x3) unit cell model 


(a) Axial Stress 


Figure 9 Stress contours for a two dimensional model of a plain weave 

composite under bending, (nominal axial strain at top surface = .00 1 ) 
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(ii) Exterior Unit Cell of a (3x3) unit cell model 


(c) Shear Stress 


Figure 9, Concluded. 
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Chapter. 5 

EVALUATION OF HOMOGENIZATION FOR GLOBAL/LOCAL STRESS 
ANALYSIS OF TEXTILE COMPOSITES 


Abstract 


Global/local analysis is essential for textile composites because of their unusually large 
microstructure. Homogenized engineering properties were used in this study to obtain global 
solutions. The response of a local region was approximated by several fundamental strain or 
stress modes. The magnitudes of these modes, which were determined from the global solutions, 
were used to scale and superpose solutions from refined analyses of the fundamental modes, thus 
obtaining a refined local solution. Results from numerical experiments showed that the use of 
homogenized engineering properties often results in significant errors in prediction of global 
response, especially at boundaries. Also, the local predictions were very sensitive to the choice 
of fundamental modes. 


Introduction 


Recently there has been an increased interest in textile composites because of potential 
increases in damage tolerance and decreased cost relative to tape laminates. These composites 
consist of a textile preform which is impregnated with resin. The interlacing used in making a 
preform can be accomplished by weaving, braiding, or knitting. Figure l shows examples of two 
weave architectures: a plain weave and a 5-hamess satin weave. (The resin pockets are removed 
in the figure so that the fiber tows can be seen.) Textile composites all have^very large 
microstructure compared to traditional tape laminates. In fact, the "microstructure can be of 

the same scale as some of the structural dimensions. 

One of the techniques proposed for analyzing textile composite structures is to use 
homogenized engineering material properties or some other measure of effective properties for 
a global analysis. This avoids the impossible burden of modeling the microstructure discretely 
in a structural model. To determine the details of the stress and strain distributions, subsequent 
analyses are performed using a refined model of a representative unit cell. The boundary 
conditions for these subsequent analyses are determined from the results of the global analysis. 
Such analyses have been discussed previously (e.g.. References 1-4). This multi-level procedure 
could be considered a global/local method and will be referred to as such herein. References 1 
and 3 discussed the accuracy of this procedure if one uses special elements (referred to as macro 
elements) for the global analysis. However, this author is not aware of any study which 
evaluated the accuracy of a global/local procedure for textile composites based on using 
homogenized engineering properties for the global analysis. 
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The objective of this paper is to describe two global/local procedures which use 
homogenized engineering material properties to expedite global stress analysis of textile 
composites and to determine the errors which are inherent in such analyses. One of the key 
questions is whether the use of homogenized engineering properties is adequate when the 
microstructure is large. To simplify the discussion and numerical experiments, only 
two-dimensional models will be examined. Admittedly, textile composites are fully 3D in their 
geometry, but the trends determined from 2D models are expected to be qualitatively correct. 

In the following sections the theoretical basis will be described first. Then the 
configurations studied will be described. Finally, the results of the numerical experiments will 
be discussed. 


Theory 

This section will describe the global/local procedures used. Figure 2 shows a schematic 
of the global/local analysis procedure. In this sketch the shading identifies the region which will 
be analyzed further using a local model. The region to be analyzed using a local model is shown 
isolated from the rest of the global model. After completing the global analysis, the boundary 
nodal displacements (u„ V;) and forces (F;, Fp are known. This boundary information is used 
to determine the appropriate loading conditions for a refined local model. There are many 
possibilities for determining these boundary conditions. In this particular study the boundary 
information was used to quantify the magnitudes of selected fundamental strain or stress modes. 
Details of the various steps are discussed in the following subsections. First, the term 
homogenized engineering properties will be defined. Then the fundamental macroscopic strain 
and str ess modes will be described, including an explanation of how the magnitude of the modes 
were determined. 

Homogenized Engineering Properties 

A unit cell is the basic building block which can be used to synthesize a woven 
composite. In this paper the woven mats are stacked symmetrically, so the unit cell consists of 
one wavelength of two mats. Homogenized engineering properties for use in the global analysis 
were determined by analyzing an infinite array of unit cells subjected to macroscopically constant 
stress states. Hence, every unit cell in the array experiences the same deformation. Periodic 
boundary conditions were applied to a single unit cell to make it behave as though it was 
embedded within an infinite array. Details about the periodic boundary conditions can be found 
in Reference 5. The homogenized engineering properties were obtained by equating energies in 
the homogenized medium to that in the actual unit cell. 

Fundamental Macroscopic Strain and Stress Modes 

In the current study the local model consisted of a refined mesh of a unit cell. In general, 
the local model could be smaller or larger. The loading for this refined unit cell was determined 
from the nodal displacements (u\ v*) or forces (Fp F‘) in the global model at the nodes which 
surround the region of interest. The local model typically has many more nooes along the 
global/local boundary than the global model. Hence, the dimensionality of the local model along 
the global/local boundary must be reduced. One technique to reduce the dimensionality is to 
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limit the response to a few macroscopic strain or stress modes. In tills paper the response of the 
local region was characterized in terms of five strain or stress modes. These modes are: 

Strain modes: 

e® : constant macroscopic c x 

e ® : constant macroscopic e y 

e® y : constant macroscopic e iy 

c® : constant gradient of macroscopic e x with respect to y 

c ® x : constant gradient of macroscopic e y with respect to x 

Stress modes: 

a° z : constant macroscopic o x 

o y : constant macroscopic o y 

o° : constant macroscopic a xy 

a® : constant gradient of macroscopic o x with respect to y 

o® x : constant gradient of macroscopic a y with respect to x 

There are interior and exterior versions of some of these modes. There are neighboring 
unit cells on all sides for interior modes and on only two sides for exterior modes. Figure 3 
shows deformed finite element meshes which illustrate the five interior stress modes. The 
shaded rectangles indicate the original mesh size and shape. The interior modes were used for 
analyzing interior cells. A mixture of interior and exterior modes were used for analyzing 
exterior cells. The mix is listed below for displacement (strain modes) and force (stress modes) 
based superposition. 


Strain Modes 

Stress 

Modes 

Mode 

Version 

Mode 

Version 

0 

C x 

interior 

0 

exterior 

0 

s 

interior 

0 

interior 

0 

interior 

0 

°xy 

interior 

0 

C *.y 

exterior 

0 

°x.r 

exterior 

0 

S.x 

interior 

0 

interior 


Only a few exterior modes were used. This is because the free surface of the exterior cell was 
a y = constant line. Some exterior modes, such as a <r y ° mode, do not exist for such a cell. 

The technique for imposing boundary conditions for the various modes is described in 
References 5 and 6. The techniques used to determine the magnitudes of the modes is discussed 
in the following two sections. 
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Strain Mode Superposition The global/local displacement field was assumed to be 
describable by the following bi-linear approximation in x and y. 


u=a+bx+cy+dxy 
v=e +fx+gy+hxy 


(1) 


The eight constants a-h can be determined by requiring that equation 1 match the displacements 
at the comer nodes of the local region. The macroscopic strain modes can be obtained by 
differentiation of the equations. The equation for e^. was further simplified by evaluating it at 
the unit cell centroid and taking it to be constant for the entire global/local boundary. This 
resulted in five strain modes: c° y , e° iy , e\ y and e° yt . In particular. 


€ x = - b *dy 

e T = g+hx 

£ v = c *f 


( 2 ) 


The coefficients b, d, g, h, and c+/are the magnitudes of the five fundamental strain modes. 


Stress Mode Superposition This technique is similar to the strain mode superposition 
method. In this case the nodal forces from the global analysis are used to determine the 
magnitudes of five fundamental stress modes. These fundamental modes were described earlier. 
This section will describe how to determine the magnitudes of these modes. 

The first step is to express the tractions T x and T y acting along the global/local boundary 
in terms of the stresses. 


T x = °x"x + °yt n j 


( 3 ) 


The relationship between these tractions and the equivalent nodal forces for a single 
element can be derived using the principle of virtual work. The result is 


F; = fT/tfiDdS 
F‘ = fcNtfdS 


where 


= 1 , number of boundary nodes 
Nj = interpolation fiinctions 


In this paper the local region is rectangular and aligned with the global xy axes so dS is 
either dx or dy. The total nodal forces for each node along the entire boundary are obtained by 
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summing the contributions from each element. Next the average stresses for the entire local 
region are assumed to be given by 


<j f * a ♦ by 
c y =c+dx 


( 5 ) 


O 


n 


- c 


These expressions for stresses are used in equations 3 and 4 to determine the equivalent 
nodal loads. Since there are many more known nodal forces (and hence more equations) than 
unknown coefficients (a-e), a least squares procedure is used to solve for the unknowns. 

Once the coefficients are determined, they are used to scale and superpose the 

fundamental stress modes described earlier. 


Configurations 


A very stubby beam was subjected to three types of loading: constant moment, distributed 
transverse shear at the end, and distributed transverse loading along the lower surface. More 
precisely, the conditions were: (see Figure 4) 


Constant moment: 

«(0,0) = v(0,0) = o 

— (^U5y) = -.01 

dy 

— (4.5,y) = -01 

dy 


Transverse end load: 

«(0jr) = v(0,y) = 0 
^(4.5 j>) = constant 


Distributed lateral load: 

u (0,y) = v(0,y) = 0 
T y (x,-3) = constant 


The beam consisted of 3x3 array of unit cells. The ratio of wavelength to mat thickness 
gives a measure of the waviness of the fiber tows. In this study this ratio (Xfh) was 1/3. 

The following material properties were assumed: 


Fiber tow (Ref. 6] 

E x = 206.900 GPa 

E = 5.171 GPa 

E x = 5.171 GPa 

v =0.25 
j y 

v = 0.25 
v* = 0.25 

zx 

G =2.386 GPa 

G =2.386 GPa 
r* 

G =2.386 GPa 

tx 


Matrix pockets 

E x = 3.45 GPa 

E y = 3.45 GPa 

E c = 3-45 GPa 

v =035 

= 035 

v" =0.35 

= 138 GPa 

G k - 1-28 GPa 

G = 1.28 GPa 
0 


Homogenized properties 
E x = 36.494 GPa 
£ y = 5.225 GPa 
E i = 36.494 GPa 
v = 1.078 
=0.154 
v K = 0.154 
G = 3.145 GPa 
G = 3.145 GPa 
G =2.000 GPa 

IX 


Figure 4 shows typical meshes which were used in this study. The reference mesh used 
5041 nodes and 1728 bi-quadratic elements to model nine unit cells. The homogenized property 
mesh used 217 nodes and 36 bi-cubic elements. The refined local mesh had 593 nodes and 192 
bi-quadratic elements. The shading indicates the two unit cells (one interior and one exterior) 
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which were analyzed using global/local analysis. Obviously, there are far fewer equations 
involved in the global/local analysis than in the conventional analysis used to obtain a reference 
solution. 


Results and Discussion 

The errors in a global/local analysis are the cumulative result of errors at the various 
stages in the procedure. To improve on a procedure requires that one know where errors are 
being introduced. Accordingly, the following discussion will begin with an evaluation of the 
predicted global response and finally examine errors in the predicted local stress distributions. 

To help evaluate the accuracy of the global analysis, the deformation of the reference 
and homogenized property meshes were compared. Figure 5 shows deformed finite element 
meshes for the three load cases. The meshes are overlaid to aid the comparison. The 
inhomogeneity in the reference mesh causes local distortions which should not (and do not) occur 
when homogenized properties are used. In Figure 5a (for a constant moment) the agreement 
appears excellent, except for the local distortion. This apparent accuracy is an artifact of the 
loading, which consisted of specified normal displacements on the left and right sides. The strain 
energy (and required moment) in the homogenized property mesh is 40% too large. In Figures 
5b and 5c the loading consisted of specified forces. The agreement between the meshes is fair 
for these cases Comparison of the strain energies in the reference and homogenized property 
models gives a scalar measure of the agreement in the predictions- The error in strain energy 
for the entire model was quite small (-6.6% for the transverse end load case and 2.6% for the 
distributed lateral load case). Also shown are magnified views of one interior and one exterior 
unit cell for each load case. (See Figure 4 for the location of the cells.) To expedite the 
comparisons, the rigid body motion of the unit cells was subtracted before plotting. Removing 
the rigid body rotation permits the unit cells to be aligned for comparison. When removing the 

du dv . 

rigid body rotation, it is important that the linear definition of rotation (l.e. rotation = — ) 

be used. For example, consider the beam in Figure 6. The beam was subjected to a moment 
at the right end. Contrary to appearances, all the unit cells have the same strain distribution. 
If the rigid body rotation is removed using the linear rotation formula, the deformed meshes for 
each unit cell will also be identical. 

The errors in the strain energies for the individual cells are tabulated below : 



Constant Moment 

Transverse End 
Load 

... 

Distributed Lateral 
Load 

Interior 

Exterior 

Interior 

Exterior 

Interior 

Exterior 

Reference 

59520 

415860 

613764 

mM 

160320 

68400 

Homogenized 

46152 

599960 

513120 

394740 

136956 

66600 

Error (%) 

-22 

44 

-!6 

-15 

-15 

-3 


The simplicity of the loading in some cases allows one to explain the source of the errors. The 
-22% error for the interior cell of a beam subjected to constant moment resulted from the 
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effective extensional modulus (which is what was used) being 22% smaller than the effective 
flexural modulus. For the interior cell of a beam subjected to transverse end load, there is both 
flexure and shear. The shear contribution to strain energy is calculated accurately but the flexure 
contribution is again low by 22%, which resulted in an net error of -16%. For the exterior cell 
of a beam subjected to constant moment the dominant deformation mode is extension. There is 
also some flexure. The 44% error in strain energy resulted from using the effective extensional 
modulus (which is based on infinite array analysis) throughout. In reality, the extensional and 
flexural modulus for exterior cells is much smaller than the effective extensional modulus for 
an infinite array. These errors illustrate the problems in using effective engineering properties 
for this class of materials. 

As discussed in the theory section, the nodal displacements and forces were used to 
determine the magnitudes of the fundamental strain and stress modes, respectively There is 
inherently some error in this approximation, regardless of the accuracy of the global analysis. 
This is because in general the actual behavior cannot be matched by just the modes selected. 
However by calculating the magnitudes of the modes using the reference mesh, one obtains a 
baseline approximation which is about as good as can be expected. Table l summarizes the 


For pure bending the strain modes for the interior cell are identical for the two meshes, 
but this is not a sign of accuracy, since the specified displacement loading required this identity. 
There was a -33 % error in the constant € y mode for the exterior cell. The other two non-zero 
modes were exact, which again was due to the boundary conditions. The error in the stress 
modes depended on the location of the cell and the particular mode. The importance of a 
particular mode cannot be seen in Table 1 . The numbers in these tables are used to scale the 

stress distributions from the fundamental solutions, i.e., o,. = c*o“ where c* = magnitudes in the 

table and = stress distribution for the "a" mode. Both the and o* must be considered 
when determining the dominant modes for a particular load case. The dominant stress mode for 
the interior cell was the gradient of a x mode, which was off by -10%. In contrast, the dominant 
mode (constant <rj for the exterior cell was off by 30%. For the transverse end^ad case the 

dominant modes were e° v , c%, and <j° v for the interior cell and 

for the exterior cell. The largest errors in the dominant modes were for t° v and o° v . These 
errors tended to be quite large. For the distributed lateral load case most of the modes were 
significant. (The c y-r and o°„ modes were not significant.) The errors in the modes tended to 
be larger than for the other two load cases. 

The magnitudes of the modes in Tables 1(a) and 1(b) can be used to scale and superpose 
displacements for the fundamental modes. These superposed displacements were determined for 
interior and exterior cells. The deformed meshes are shown m Figures 7 and 8 for strain and 
stress mode superposition, respectively. As was done in Figure 5, the rigid body components 
were subtracted to make the comparisons of the unit cell deformations more accurate. The 
thicker lines indicate the superposition results. The results labeled "Reference Superposition 
were obtained by using the reference mesh to determine the magnitudes of the modes 1 he 
"Reference Superposition" results show that even if a global analysis is exact, the loca 
deformation cannot in general be represented in terms of a few fundamental modes. Regardless 
of the type of loading, the interior behavior is more closely approximated than the exterior 
behavior. Strain mode superposition appears to be more accurate for interior cells. In contrast. 
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stress inode superposition appears to be more accurate for exterior cells. 

The next step was to determine the accuracy of the calculated stresses. Three types of 

solutions were examined: 

1) the reference solution 

2) the global/local solution in which the global analysis used the reference mesh to 
determine the modal components 

3) the global/local solution in which the global mesh used homogenized material 
properties. 

The results are presented two ways. First, just the peak stresses for each load case and analysis 
type will be summarized in tabular form and then a few stress contour plots will be discussed. 

The peak stresses are tabulated in Table 2. The errors in the global/local stress 
calculations varied over a wide range. The simplicity of the loading for the constant moment 
case eliminated one source of error... that related to determining the average strain field using 
five strain modes. The second potential source of error was in determining the local stresses 
from the fundamental strain modes. This was no problem for the interior cell; the error was 
essentially zero. In Reference 7 it was shown that unit cells at least one cell away from a free 
surface behaved very much like ones embedded in an infinite array. Since the fundamental strain 
modes (which included two bending modes) were based on infinite arrays, the accurate 
prediction is no surprise. In contrast, the errors are significant for the exterior cell. Even when 
the refined reference mesh was used to determine the magnitudes of the different modes, the 
errors were not negligible. The stress mode superposition method tended to perform better for 
exterior cells than the strain mode method. The response of an exterior cell is complex and 
hence poorly represented by the particular few strain or stress modes considered. The errors due 
to modal reduction increased with the complexity of the applied load. For the distributed lateral 
load case the errors were significant for the interior cell and intolerable for the exterior cell. 

Obviously, the interior and exterior unit cells experience different loading and different 
modes are dominant for the two cells. Hence, the larger errors for the exterior cell could be 
due to errors associated with particular modes, rather than the location of the cells. This was 
checked in an approximate sense by adding a layer of unit cells to the top and bottom of the 
current global model. Global/local analysis of this thicker beam was performed for transverse 
end load case. The errors in the peak stresses for the unit cell which had been on the exterior 
for the thinner model were now much less. This suggests that the behavior of an exterior cell 
is inherently more complicated than that of an interior cell. 

Examination of errors in predicted peak stresses gives only a limited appreciation of the 
accuracy (or inaccuracy) of the predictions. Figures 9 and 10 show stress contours for the 
transverse end load case for interior and exterior cells. The contours for the interior cell (Figure 
9) for the global/local analysis match very closely with the reference solution. The contours for 
the exterior cell (Figure 10) are not as close, but still seem to agree fairly well, even though the 
errors in the peak stresses are up to 26%. Peak stresses will probably not be useful for 
predicting failure, since they occur at a point (or at least a very small region). A critical stress 
criterion will probably have to consider the average stress in some characteristic volume. The 
visual similarity of the contours in Figures 9 and 10 suggests that when global/local analysis is 
used, the errors in a practical failure criterion might not be as bad as the errors in a peak stress 
criterion. This visual evaluation of the similarities in the contours in Figures 9 and 10 is 
subjective and could be wrong. A more objective method is badly needed. 

One technique which was considered was plotting the tow area which had a stress greater 
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than a particular value. The reasoning is that if a large stress occurs over only an extremely 
small region, then the stress calculation is suspect, since the scale is too close to that of the fiber 
diameter. Figures 1 1 and 12 show results for a, for an interior and exterior cell, respectively 
for the transverse end load case. Two graphs are shown in each figure : one which includes the 
entire stress range, and the other which zooms in on the peak stress region. Results from the 
four analyses agree quite well for the interior cell. Although the homogenized superposition 
technique has a 15 percent error in predicting the peak stress, the distribution is predicted fairly 
well If the failure criteria requires that a particular volume be at a critical stress, the prediction 
would be in error about 15 percent for very small critical volumes, but the error would be less 
for larger critical volumes. Figure 12 shows analogous results for the exterior cell. The errors 
are larger than for the interior cell, but the trend is the same, ie. for larger critical volumes the 
error in failure prediction is less than for very small critical volumes. 

Concluding Remarks 

Global/local stress analysis techniques based on the use of homogenized properties for 
the global analysis were evaluated. A very stubby beam containing nine unit cells was subjected 
to three types of loading. Considering the strong macroscopic stress and strain gradients relative 
to the microstructure these were probably fairly severe tests. For force type loading the overall 
deformation of the beam was not always predicted very well using homogenized properties. For 
larger configurations with more unit cells (and hence more homogeneous microstructure) the 
accuracy is expected to be considerably better. The accuracy of the calculated stresses was not 
too bad for interior cells, but was poor for exterior cells. This is not surprising based on earlier 

work on free boundary effects. _ J .._ w . 

Regardless of how a global solution is obtained, there is considerable difficulty in using 
the crude nodal force and displacement information from the global mesh to determine 
appropriate load conditions for the local mesh. In this paper a modal technique was used. For 
the constant moment and transverse end load cases this technique performed well. For the more 
complicated case of distributed lateral load the performance was only fair for the interior cell 
and poor for the exterior cell, even when a refined global mesh was used. 

There are several steps (and inherent approximation at each step) in global/local analysis. 
This study was just a beginning. Further work is needed in several areas. Alternatives to 
homogenization , such as the macro elements in References 8 and 9 need to be evaluated. Other 
techniques for imposing the global solution on a local model also need evaluabon, including 
additional types of fundamental modes and the use of smaller local regions. Finally, more 
realistic configurations need to be identified and studied. Otherwise it is difficult to assess the 
significance of the errors in the various global/local techniques for practical applicauons. 
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Table 2 Reference peak stresses and corresponding global/local stresses (». - P«ceni .fro.) 
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Figure 1 Examples of textile architecture. 



82 




83 


3 Deformed meshes for the five interior stress modes. 






9.0 


(a) Reference mesh 



(b) Mesh used with homogenized properties 


X 



(c) Local finite element mesh 


Figure 4 Finite element meshes. 
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(a) Constant Moment 


(b) Transverse End Load 


(c) Distributed Lateral Load 
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Interior Cell 



Interior Cell 



Interior Cell 


Figure 5 Comparison of deformed meshes of reference model and 
homogenized property model. 
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Figure 6 Artifacts due to magnification of deformation. 
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Homogenenized Superposition 



Exterior 


(a) Constant Moment 



Interior 


(b) Transverse End Load 




Exterior 





Figure 7 Comparison of deformed meshes for strain mode superposition results. 
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Figure 8 Comparison of deformed meshes for stress mode superposition results . 
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Reference solution 


Homogenized superposition 




+ 1300c +08 
+l.200e+08 
+9.000c+07 
+6 ,000c +07 
+3.000e+07 
+0.000e+00 
-3,000c+07 
-6.000c +07 
-9.000c +07 
-1.200c+08 
-l.SOOc+08 



Peak magnitude : 2.47E+08 


Peak magnitude: 2.I1E+08 

Error : -15% 



Peak magnitude: 6.10E+07 


a 


y 



+4.000e+07 
+3.200e+07 
+2.400C+07 
+ 1.600c +07 
+8.000e+06 
+0.000e+00 

-8.000C+06 

-1.600c+07 

-2.400e+07 

-3.200c+07 

-4.000e+07 



Peak magnitude : 5.13E+07 

Error : -16% 



Peak magnitude: -1.15E+08 


Oxy 

+8.750e+07 
, , +7.000e+07 

[gj — +5.250e+07 
m +3300e+07 
+1.750C+07 
+0.000e+00 
-l.750e+07 
-3.500e+07 
-5.250C+07 
|- -7.000c+07 
-8.750e+07 



Peak magnitude : -0.99E+08 

Error -14% 


Figure 9 


Comparison of stress distributions of reference solution and homogenized strain mode 
superposition for the interior cell for transverse end load case . 
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Reference Solution 
Reference Superposition 
Homogenized Superposition 
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Figure II Tow area with a x greater than particular stress level 
Interior cell for transverse end load case. 



Reference Solution 
Reference Superposition 
Homogenized Superposition 



(%) 29JV ^01 
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Figure 12 Tow area with Ox greater than particular stress level 

Exterior cell for transverse end load case. 



Chapter: 6 

Simulation of Progressive Failure in Plain Weave 
Textile Composites 


Abstract 

Three-dimensional finite element analysis was used to simulate progressive failure 
of a plain weave subjected to in-plane extension. The loading was parallel to one of the 
tow directions. The effects of various characteristics of the finite element model on 
predicted behavior were examined. More numerical studies and comparisons with 
experimental data are needed to establish guidelines for accurate progressive failure 
prediction. Also the sensitivity of the predictions to the tow waviness was studied. The 
predicted strength decreased considerably with increased waviness. 

Introduction 

Textile composites consist of interlaced fiber bundles which are then impregnated 
with a matrix material and cured. Figure 1 illustrates the architecture for a plain weave 
composite. The interlacing of the fibers offers the potential for increased through- 
thickness strength There is also the potential for reduced fabrication costs, since fairly 
complicated shapes can be formed using textile machinery. One disadvantage of textiles 
is the difficulty in predicting their performance. The complex geometry makes detailed 
stress analysis quite challenging. The early analyses were based on modified laminate 
theory, (eg. References 1,2). In recent years there have been a few attempts to 
discretely model the fiber bundle architecture and predict internal stress states (eg 
References 3-10). Reference 10 presented a particularly interesting progressive failure 
analysis of a plain weave composite. The results in Reference 10 consisted of nominal 
stress strain curves. The response of the composite was almost linear for in-plane 
extension and highly nonlinear for in-plane shear, the nonlinearity was primarily a result 
of progressive damage. However, little information was provided on damage evolution 
and load redistribution within the composite during the loading process. Also, there was 
no indication of the sensitivity of the predictions to mesh refinement or other approxima- 
tions inherent in such analyses. 

This paper has two objectives. The first is to evaluate the sensitivity of predicted 
progressive failure to quadrature order, mesh refinement, and choice of material 
degradation model. The second objective is to describe the nature of the progressive 
failure process for two weaves with very different waviness. Loading consisted of a 
nominally uniaxial stress along one of the fiber tow directions. Only mechanical loads 
were considered in this study. To simplify the response the composite was assumed to 
consist of an infinite number of unit cells in all three coordinate directions. 
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The following sections begin with a description of the basic theory used for 
progressive damage modelling. Then the configurations will be described. Finally the 
results from the numerical simulations will be discussed. 

Theory 


There is no "right" way to model damage evolution that is also practical. It is not 
feasible to discretely model the damage, so approximation is unavoidable. Perhaps the 
simplest procedure to account for damage in a finite element model is to modify the 
constitutive matrix at the quadrature points of a numerically integrated finite element. 

No history effects are included, so the analysis of the loading becomes a series of elastic 
analyses. Of course, there are many possibilities for how to modify the constitutive 
matrix. Three techniques were used herein. The first method considered the material 
totally failed (i.e. the entire constitutive matrix was reduced to essentially zero) when any 
allowable stress component was exceeded. This method will be referred to as the non- 
selective discount method. Except as noted, this technique was used in the analyses. 

The second technique selectively reduced the rows and columns of the constitutive matrix 
according to the particular stress allowable which was exceeded. The third technique 
selectively reduced the engineering moduli according to the particular stress allowable 
which was exceeded. The scheme for this selective reduction was based on References 
10 . 


Figure 2 gives a flowchart for the progressive failure analysis. First a linear 
analysis was performed. Based on the calculated stresses, the initial load was scaled back 
so that failure would occur only at points which were within two percent of the maximum 
normalized stress. (The stresses were normalized by the respective strengths.) The 
constitutive matrix was modified at the failure points. Residual forces were calculated 
and used to determine the incremental displacements required to restore equilibrium. 
The total displacements were updated and used to determine the new stresses. If no 
further failures occurred at the current nominal strain state, the nominal strain was 
incremented to cause failure. This procedure was repeated until there was total failure 
or at least loss of most of the original stiffness. 

Configurations 

The fiber bundles or tows in the models were generated by translating a lenticular 
cross-section along a sinusoidal path. The waviness ratio is defined to be the ratio of the 
woven mat thickness to the wavelength. Except where indicated otherwise, the results 
presented are for a waviness ratio of 1/3. More details about the mesh geometry can be 
found in Reference 8. The following subsections describe the finite element meshes, the 
boundary conditions, and the material properties. 
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Finite Element Meshes 


Symmetry in the material and loading was exploited so that only 1/32 of a unit cell 
had to be modeled. A wide range of mesh refinements were used, as shown in Figure 3. 
The crude mesh had only 4 elements and 43 nodes. The most refined mesh had 192 
elements and 1049 nodes. 


Boundary Conditions 

The periodic boundary conditions for a complete unit cell are quite simple. The 
appropriate boundary conditions for a 1/32 unit cell are a bit more complicated. Deriva- 
tion of the periodic boundary conditions is somewhat tedious, so details will not be given 
here. Details can be found in Reference 8. The periodic conditions are listed below. 
Figure 3 shows the coordinate system assumed. 


u(a/2,y.z) = Uq 
u(0,y,z) = -u(0,y.-z) 

u(x,0,z) = u(x,0,-z) 

w(x.y.c/2) as constant 
w(O.y.z) = -w(O.y.-z) 
w(x.O.z) = -w(x.0.-z) 


v(x^/2,z) = constant 
v(0,y.z) = v(0.y,-z) 

v(x.0.z) = -v(x.0.-z) 


The load was controlled by specifying the magnitude of i^. 


Material Properties 

The unit cells contains two "types of materials: the tows and the matrix pockets. 
Relative to the material coordinate system, the properties of the tows are invariant 
(before damage occurs). Of course, the properties of the tows are needed in the global 
coordinate system. Fourth order tensor transformation formula were used to perform 
the required calculations. The rotation angles to be used in these formulas were 
obtained at each quadrature point by using interpolation. This procedure was shown in 
Reference 8 and 11 to be preferable to using a single angle for the entire element. The 
particular properties used are listed below. These properties are from Reference 12. 
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E.i 

E22 

E33 

G l2 

0,3 

G 2 3 

V l2 

Vl3 

V23 


Tow properties 


Modulus 

Strength 

154.27 

GPa 

2342.0 

MPa 

10.80 

GPa 

56.6 

MPa 

10.80 

GPa 

56.6 

MPa 

7.47 

GPa 

48.7 

MPa 

7.47 

GPa 

48.7 

MPa 

333 

GPa 

48.7 

MPa 


0.278 

0.278 

0340 


Matrix properties 


Modulus 

Strength 

3.45 GPa 

84.85 MPa 

3.45 GPa 

8435 MPa 

3.45 GPa 

84.85 MPa 

138 GPa 

101.00 MPa 

138 GPa 

101.00 MPa 

138 GPa 

101.00 MPa 

035 


035 


035 



Results and Discussions 

Most of the results in this paper illustrate the effects of characteristic of the finite 
element model on the progressive failure prediction. The effects of quadrature order, 
mesh refinement, and material degradation strategy will be considered first. Then the 
effect of tow waviness on failure behavior will be discussed. 


Figure 4 shows the effect of quadrature order on the stress-strain curve. The 
peak stress obtained using 8 quadrature points (2x2x2), is 10 percent higher than that 
obtained using 27 or 64 points. Although the peak stress is the same for 27 and 64 
points, damage is predicted earlier when 64 points integration is used. This sensitivity is 
not particularly surprising for at least two reasons. First, when more quadrature points 
are used, the more extensive sampling is more likely to find the extremes in the stress 
field. Second, when failure occurs within an element and the constitutive matrix 
is modified, the element becomes inhomogeneous. The numerical integration effectively 
fits a polynomial function to the variation of material properties. Since the properties 
are very different in the failed and unfailed parts of the element, it is difficult to obtain a 
good fit. In fact, there is concern as to whether the assumed quadratic displacement 
functions for a 20-node element are sufficient to obtain a reasonable approximation 
regardless of the integration order. 

Figure 5 shows the effect of mesh refinement on the predicted stress-strain curve 
for two waviness ratios. The 4 element model predicts the correct trends, but is quite 
inaccurate. The error is much worse for the large waviness ratio. For the 1/6 waviness 
ratio, the 32 and 192 element models agree quite well. There is a considerable differ- 
ence between the 32 and 192 element models for the 1/3 waviness ratio. Although the 
response is quite brittle for both waviness ratios, there is no more non-catastrophic 
damage before collapse for the larger waviness ratio. 
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Figure 6 shows the effect of the discount factor on the stress-strain curve. The 
stiffness terms at failed quadrature points were reduced to either 01. or .0001 of the 
original value. Intuitively, one might expect to obtain the same results. Figure 6 shows 
that there was no difference in the peak stress, but the response was very different when 
there is considerable damage. 

Figure 7 shows the stress-strain curves obtained using non-selective discount 
method, selective reduction of rows and columns in the constitutive matrix (the stiffness 
terms, not the compliance terms), and the selective method described in Reference 10. 
The selective method described in Reference 10 predicts about a 21 percent higher peak 
stress than the non-selective discount method. 

Figures 8 and 9 show the effect of mesh refinement and waviness ratio on damage 
accumulation during loading. The black region indicates the damage zone. The stress- 
strain curve for a particular mesh is shown above the results for that mesh. The points 
labeled A,B, and C, indicate the correspondence between the strain level and the damage 
contours. Also indicated are the stress components which contributed to the damage 
contour. The 4 element mesh does not perform well at all for the waviness ratio of 1/3, 
but does a little better for the 1/16 waviness. The 32 element mesh performs reasonably 
well for obtaining qualitative results. Further numerical studies are needed to determine 
how close the 192 element mesh is to convergence. For the 1/3 waviness ratio the a 33 
stress component dominates the damage development up to the point shown. For the 
1/6 waviness ratio cr 33 plays a part, but there is also significant cracking of the 90 degree 
tow due to 022 • Reference 6 had also noted a change in initial damage mode with 
waviness ratio. 


Concluding Remarks 

Simulation of progressive failure in a plain weave composite is extremely complex. 
Consequently, only approximate treatment is practical at this time. One of the goals of 
this paper was to examine the effect of several approximations on predicted behavior. 
The one obvious conclusions from this study is that the predictions are quite sensitive to 
a number of decisions which must be made when assembling a finite element model. 
Further numerical experiments and comparisons with experimental data are needed to 
establish guidelines for accurate analysis of progressive failure. Another objective of this 
paper was to describe the effect of tow waviness on damage accumulation. The results 
suggest that the degree of waviness not only affects the stress at which damage initiates, 
but also the type of damage which occurs. 
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Figure 2 Flowchart of progressive failure analysis. 
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(a) Waviness ratio = 1/3 

Figure 5 Effect of mesh refinement on stress-strain response 
Non-selective discount method was used. 
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Figure 7 Effect of material degradation method on stress-strain response. 













Chapter: 7 

Effect of Assumed Tow Architecture on Predicted Moduli 
and Stresses in Plain Weave Composites 


Abstract 

This paper examines the effect of assumed tow architecture on the predicted 
moduli and stresses in plain weave textile composites. Two architectures are examined 
which have a sinusoidal tow path and a lenticular cross-section. Three-dimensional finite 
elements are employed to model a T3 00/Epoxy plain weave composite with 
symmetrically stacked mats. Macroscopically homogeneous in-plane extension and shear 
and transverse shear loadings were considered. Symmetries are exploited which 
permitted modeling of only l/32nd of the unit cell. Accounting for the variation of 
material properties throughout each element is determined to be necessary for accurate 
prediction of stresses in the composite. For low waviness, the two tow architectures 
examined are very similar. At high waviness, the stress predictions are much more 
sensitive to the assumed tow geometry. 


Introduction 

As more demands are placed on structural materials, more complex materials 
must be developed in order to satisfy these demands. Fibrous composites are being 
employed for many of these applications. One type of fibrous composite which is 
receiving increased attention is the woven compodte, which is constructed by grouping 
the fibers into bundles called tows and weaving these bandies together to make preforms. 
With the addition of a matrix material and curing, a woven composite is constructed. 
Figure 1 shows a schematic of a plain weave composite. 

Woven composites have recently been receiving considerable attention due to their 
potential for improved properties compared with traditional laminated composites. Some 
of these properties include an increase in impact toughness and higher specific-strength 
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and specific-stiffness [1]. Other properties include the ability of near-net-shape 
manufacturing capability [2,3] and a high resistance to damage and impact [4]. 

Unfortunately, the microstructure of a woven composite is highly complex making 
it usually not cost effective to experimentally explore the effect of different weave 
geometries [1]. In fact, it is not currently practical to model the actual tow geometry. 
Idealization is mandatory. The idealization is critical, since die composite properties are 
very sensitive to some aspects of the tow geometry [5]. Since the stress analysis problem 
is too complicated for the luxury of modeling insignificant details, the obvious goal is to 
include just enough details in the idealized geometry to predict die most important aspects 
of the composite behavior. One way to determine the necessity of modeling particular 
details is to perform analytical studies which use various approximations of the tow 
geometry and observe the effects of these approximations on predicted response (eg. 
moduli and internal stresses). These analytical studies should be supplemented by 
experiments. However, this paper will focus only on the analytical task. 

The objective of this paper is to determine the sensitivity of predicted moduli and 
internal stress distributions and failure initiation to one aspect of tow geometry: the 
variation of tow cross-section along the length of the tow. There are basically two 
categories of stress analysis methods which could be used for the study. One of these 
is based on a modified form of laminate theory (eg. References [4,6-11]) and the other 
is finite elements (eg. References [1,5,12-18]). The first category of analyses is only 
suitable for modulus prediction and rough estimates of stresses. To fully evaluate the 
accuracy of various approximations of tow geometry requires numerical analysis. Finite 
elements are ideal for this task and were used in this study. 

In this study, only elastic behavior of mechanically loaded specimens will be 
considered. Residual stresses caused by differences in coefficients of thermal expansion 
should be included in a later study. 

Although there are many types of woven composites, plain weaves will be studied 
in this paper because it offers sufficient complexity for the task at hand. The tows of a 
plain weave are woven together as shown in Figure 1. The tows which run in the 
x -direction are the warp tows. The tows running perpendicular to these are the fill tows. 
In a balanced weave the warp and fill tows have the same geometry and material 
properties. A balanced weave was used in this study. 

The next section describes the assumed tow architectures. Next, boundary 
conditions used with the finite element model of the infinite periodic arrays will be 
described followed by a description of the material properties. Finally, the results of this 
study will be presented followed by conclusions. 
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d etermin ed for in-plane extension and shear and transverse shear. Included in this 
examination was a study on the specification of material properties using the single angle 
and multi-angle methods described earlier. 

A convergence study was performed for a WR=l/3 to determine a suitable mesh 
for prediction of elastic properties and stress distributions for the two architectures. 
Twenty node isoparametric hexahedral dements were used. The meshes used are shown 
in Figure 5 for the translated and extruded models. Convergence of the elastic properties 
for t ranslated tows is shown in Figure 6. The translated moduli are nor maliz e d with 
respect to the converged extruded modd in order to also show the relative difference in 
moduli between the two architectures. This figure also compares the multi-angle and 
single angle methods of specifying material properties in the elements. As can be seen 
from the figure, convergence of most of the elastic properties was reached with just 32 
elements using the multi-angle method. 

Convergence of the stress distributions is shown in Figure 7 for translated tows 
using the multi-angle method. This figure shows the an stress distribution in the tows 
after subjecting the meshes to a 1 % strain. Even using the 454 modd with 400 
elements convergence was not reached. However, there was only a moderate change 
between the 343 model with 192 elements and the 400 element 454 modd. For this 
reason, the 454 model will be assumed to be nearly converged and will be used to study 
both the elastic moduli and stress distributions for the translated and extruded models. 

A further comparison of the single angle and multi-angle method is shown in 
Figure 8. This figure shows the a u stress distributions in translated tows using single 
and multiple angles in the dements for a mesh subjected toa 1% « n strain. Due to 
symmetry, the stress at the end of the tow as indicated by A should be zero, but as 
shown in the figure, the predicted stress is not zero for the angle angle method. In fact, 
the highest stress in the tow using the angle angle method is reported to be at A. If the 
mesh were refined, the incorrect stresses would disappear; however, this would be 
extremely costly using 3D finite elements. This indicates that the multi-angle elements 
are essential to an accurate analysis of the composite. Therefore only the multi-angle 
method will be used in for the following analyses. 

The tow volume fraction vs. waviness ratio of translated and extruded tows is 
shown in Figure 9. (Note that this is not the fiber volume fraction.) The difference 
between the translated and extruded models increases with waviness. At WR=l/3, the 
percent difference between the tow volume fraction is approximately ,4.5%. This is 
significant because this alone will cause a difference in the predicted moduli of the 
composite. Additionally, as the waviness approaches zero, it is seen that the volume 
fraction for translated and extruded tows converge. Something else interesting to note 
is that the volume fraction of the translated tow is constant with respect to changes in 
WR_ which seems unrealistic. The reason for this is that the volume of a wavelength of 
the translated tow is simply the cross-section area times the wavelength. 
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Figure 10 shows the differences in predicted elastic properties for weaves with 
translated and extruded tows. At low waviness, the figure shows that the differences 
between predicted elastic properties for translated and extruded tows are very small. 
However, at higher waviness, the difference between extruded and translated tows is 
more apparent. The difference is greatest for E*. The translated model predicts a 6.8% 
higher moduli than the extruded model at a WR= 1/3. For WR< 1/6, the difference in 
predicted moduli is less than 1 % for all of the moduli. 

Stress distributions are not as easily compared as the moduli due to the amount 
of data which must be presented. For this reason only stresses which contribute 
significantly to the predicted failure of the material will be presented. These stresses are 
determined by normalizing the stress in a material with respect to its material strength. 
Also, since the greatest difference occurs at a WR=l/3 for the three WR’s studied, only 
stress distributions at this WR will be presented. 

Three macroscopically homogenous single stress states were considered: <r„, o^, 
and a u . These three loadings were imposed on the l/32nd sub-cells shown in Figure 5 
with the boundary conditions shown in Table 1. Stresses for the l/8th sub-cell were then 
generated using the appropriate factors given in Table 2. 

Extension in warp direction (o a loading) 

A extensional load was applied to the 454 translated and extruded meshes in 
Figure 5 using the boundary conditions shown in Table 1 so that a 1 % e a macroscopic 
strain was applied to the translated model. The same force need to apply the 1% 
macroscopic strain was then applied to the extruded model. Since the stiffness is 
different for the two architectures, the strains are different Figure 11 shows the 
predicted Tsai-Hill failure criterion for the warp tows. Only die warp tows are presented 
because failure is predicted to initiate in the warp tows. Although the contours appear 
very similar, the maximum failure criteria predicted by die extruded model is 22% higher 
than that predicted by the translated model. This indicates that the translated model can 
withstand 10% more load than the extruded before initial failure occurs. By examining 
the maximum normalized stresses experienced by the warp tows, the stresses in the warp 
tow which are the main contributors to failure were determined to be the <r n and a 33 
stresses. Consequently, only these two stresses will be examined in the remainder of this 
sub-section. 

Figure 1 1 also compares the normalized a n and a }3 stress distributions predicted 
in the warp tows with the translated and extruded models. As can be seen in Figure 11, 
cr n stress distributions in the regions of failure in the warp tows are very similar and 
contribute little to the initial failure for the translated and extruded models. However, 
near the edge of the tows one can see a more significant difference in the two models. 
The translated model predicts a high <r n stress in the cross-over region as shown in 
Figure 2. This is not seen in the extruded model in which parallel tows are separated by 
a region of resin. This difference is interesting but somewhat insignificant since failure 
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does not iniatc in this region. 

The Ojj stress distribution indicates that initial failure is being caused primarily 
by this stress. The tows are failing in the region where die warp tows in adjacent mats 
are closest to one another because the warp tows are wanting to straighten and pull apart 
from one another. This prediction should change drastically if the mats are amply 
stacked instead of symmetrically stacked as assumed in this study. The extruded mesh 
experiences a maximum Ojj stress 10.3% greater than the translated mesh. Other than 
this, there is no substantial difference in the warp tow stress distribution. 

As stated earlier, the stress levels in the fill tows and resin regions of the 
composite are uniteresting due to the relative size of the stress exhibited in these regions. 
However, if thermal stresses were included in the analysis, these regions could become 
very interesting. 

In-plane shear (a^ loading) 

A Ojy shear load was applied so that al%^ macroscopic strain was experienced 
by the 454 translated model in Figure 5. This load was then applied to the 454 extruded 
model. Graphs of the Tsai-Hill failure criteria are shown in Figure 12 for the tows. 
Note that only one tow is shown because the stress state in the warp and fill tows is the 
same. The maximum predicted failure criterion for the translated model is 30% larger 
than the maximum predicted by the extruded model. This is a substantial difference. 
Using (14), the translated model will initiate failure at a macroscopic stress level 14% 
lower than the extruded model. The normalized stresses indicate that the primary 
contributing stresses to the initial failure in the warp tow of the composite are the 
shearing stresses, <r 12 and <%. For the fill tow, the stresses are <r, 2 and <r 13 . 

The normalized a, 2 stress distributions in the tows are also shown in Figure 12 
for the translated and extruded models. The maximum stress predicted by the translated 
model is 8% greater than that predicted with the extruded model. This indicates that 
parallel tows in the composite are influencing the stresses in this region. Due to the 
extruded tows being separated by regions of resin in the cross-over region as shown in 
Figure 2, the extruded model predicts a higher strain to initial failure than that predicted 
by the translated model. 

The comparison of the normalized stress distributions in the warp tows is also 
shown in Figure 12. This figure shows that there is a considerable difference in the 
predicted stresses of the two architecturess. The translated architecture predicts 
approximately 50% greater stress than the extruded at the cross-over region in the l/8th 
sub-cell. The extruded architecture predicts its highest stress in a region away from the 
center of the l/8th sub-cell as shown in the figure. The reason for this is, again, that the 
tows of the extruded model are seperated by regions of resin. The o n stress distribution 
is what causes the large difference in the predicted failure of the two architectures. 
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Trasverse shear (cr a loading) 

A o a shear load was applied so that a 1 % £n macroscopic strain was applied to 
to the 454 translated model in Figure 5. This lad was then applied to the 454 extruded 
model. The Tsai-HUl criteria is plotted in Figure 13 for these two models. Failure is 
predicted to initiate in the resin between the mats as indicated by A in the figure. The 
translated model predicts a failure criterion 1.37% less than that predicted by the 
extruded model which indicates that initial failure will occur at essentially the same 
macroscopic stress level in both models. The main source of initial failure is the<r 13 
stress. If the mats were assumed to be simply stacked, a much differenent prediction 
would be made for this load case. 


Conclusions 

The effect of assumed tow architecture, specifically the effect of how a cross- 
section is assumed to sweep out the volume of a tow, was studied for a plain weave 
composite with symmetrically stacked mats. Comparison of two architectures, translated 
and extruded, was made on the basis of predicted moduli and stress distributions. 
Results of the analysis show that at a high waviness ratio there is a significant difference 
between the two architectures. 

It was determined early in this study that to provide an accurate analysis, the 
variation of material properties through the elements must be modelled accurately. This 
was shown by the convergence of the elastic properties and stress distributions of the 
multi-angle method compared to the single angle method. A more refined analysis was 
also achieved by reducing the size of the model which needed to be analyzed. It was 
determined that only l/32nd of the unit cell for symmetrically stacked mate need be 
analyzed. Very dramatic computational savings were obtained by this reduction. 

The differences in the tow architectures revealed several drawbacks and 
advantages of the translated and extruded models. The translated model thins at high 
waviness ratio which is unrealistic; however, the mesh generation of the translated model 
is algebraic. The extruded model does not thin at high waviness but at the cost of 
complicating the mesh generation. The extruded model also includes thin matrix layers 
between the warp and fill tows which the translated does not. The significance of this 
was not determined. 

Of the three waviness ratios studied, the largest difference was observed at the 
highest waviness, WR=l/3. Very little difference between the translated and extruded 
architectures was observed at the other two waviness ratios. Two of the loading cases, 
extension in the warp direction and in-plane shear, indicate that there is a significant 
difference between the two architectures. Extension in the warp direction indicated a 
significant difference in the predicted extensional moduli, E„ and in the load to initial 
failure. In-plane shearing showed a significant difference in the stress state in the cross- 
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over region. The third loading case, transverse shearing, showed little difference in the 
predicted moduli and stress distributions due to how the mats are assumed to be stacked. 
If the mats were assumed to be simply stacked instead of symmetrically stacked more 
difference might be noted in transverse shearing as well as extensions! loading. 
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Table 1: Boundary conditions for l/32nd sub-cell for extension, transverse shear and in- 
plane shear loadings. 


I Macroscopic extension: 


u(^^)=«o =v o 

u(Oj'^)=-u(Oo',-z) v<0,y,z)=v(0,y,-z) 

u(x,0,z)=k(x,0,-z) v(x,0,z)=-v(x,0,-z) 

"'(Oo'.z) = - h<0o', -z) 
w<x,0,z) = - h<x j O, -z) 

Transverse shear (<r u loading): 

<xj,±^)=±u 0 v(x,-|^)=0 

II 

r 

«(0,y r s)=-u(0,y,-z) v(0,y,z)=v(0 ,y,-z) 

'KQo^)=-w(0o'>-z) 

u(x,0,z) =u(x,0, -z) v(x,0,z) =v(x,0, -z) 

h<x,0,z) = -> t<x,0, -z) 

In-plane shear (o^ loading): | 

«(x,-|^)=d v(^,y,z)=d wO|,y,z)= 

=*{x,ij>e) =w(xo',±-|) =0 
2 2 

«(0,y^)=«(0,y,-z) v(0^)=vCOo',-z) 

wfOovz) = -w(Oo'» ~z) 

u(x,0^)=u(x,0,-z) v(x,0,z)=v(x,0,-z) 

h<x,0^)=-h<x,0,-z) 
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Table 2: Factors for synthesis of l/8th sub-cell stresses from l/32nd sub-cell. 
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Equations for formation of stresses in quadrants 2,3, and 4 in the xy plane. 


o?(-xy,-z) = 




o*(x,-y,-z) = o](x,y^)f* 









































































































Table 3: Comparison of CPU and memory requirements for three models with same 
refinement Times were determined on an IBM RISC System/6000 POWERstation 355. 


Cell Size 

CPU Time 
(min):(sec) 

Elements 

Memory 
Requirements 
(Words / HloBytes) 

Full Unit Cell 

37:05 

512 

7.932K / 63,456K 

S l/8th Unit Cell 

0:28 

64 

310K / 2,480K 

1 l/32nd Unit CeU 

0:10 

16 

38K / 304K 
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used to define cross-section also shown. 




Figure 2: Comparison of translated and extruded architectures. 
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Figure 4: Assumed geometry and unit cell for symmetrically stacked mats. 
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Figure 5: Meshes used in study of translated and extruded tows. 
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(c) Poisson's Ratios 

Figure 6: Convergence of predicted moduli for translated tows normalized with 
respect to extruded 454 model. 
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Figure 9: Variation of tow volume fraction with waviness ratio for extruded and translated models 




a) Extensional Moduli 




b) Shear Moduli 



Waviness Ratio 



c) Poisson's Ratios 

Translated 

Extruded 

Figu re 10: Comparison of predicted elastic moduli using translated and extruded 
models. Values are normalized with respect to values predicted at 
WR=0.05. 
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Translated 


Extruded 


Tsai-Hill 



+1375e+01 

+1^50e401 

+1.125e+01 

+1.000e-»01 

+8.750e+00 

+7500e400 

+6.250C400 

+5.000e400 

+3.750e+00 

+2.500e+00 

+1.250e+00 



<>11 



+2.354c+00 

+2.208e+00 

+2.062e+00 

+1.917e+00 

+1.771e+00 

+l.625c+00 

+1.479e+00 

+1333c+00 

+1.188e+00 

+1.042e+00 

+8.958e-01 



°33 



+4J83e+O0 

-f4.167e+00 

+3.750e+00 

+3.333c+00 

+2.917e+00 

+2.500e+00 

+2.083e+00 

+1.667e+00 

+1.250e+00 

+8.333e-01 

-t-4.167e-01 



Figure 11 : Comparison of predicted Tsai-Hill failure criterion, normalized Cj and 
normalized O33 in warp tows of translated and extruded meshes subjected 
to load generated by 1% e xx extension of translated model. 
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Translated 


Extruded 


Tsai -Hill 



+2.800e-01 

+2.600e-01 

+2.400e-01 

+2.200e-01 

+2.000e-01 

+1.800e-01 

+1.600e-01 

+1.400e-01 

+1.200e-01 

+1.000e-01 

+8.000e-02 



a 12 



+4.417e-01 

44.233e-01 

44.050e-01 

+3.867e-01 

+3.683e-01 

+3.500e-01 

+3.317e-01 

+3.133e-01 

+2.950e-01 

+2.767e-01 

+2.583e-01 


°23 



+1.833e-01 

+1.667e-01 

+1.500e-01 

+1.333e-01 

+1.167e-01 

+1.000e-01 

+8.333e-02 

+6.667e-02 

+5.000e-02 

+3.333e-02 

+1.667c-02 



Figure 12: Comparison of predicted Tsai-Hill failure criterion, normalized crj-- and 
normalized 0^3 in warp tows of translated and extruded rr.eshes sur;ected 
to load generated by 1% e xy shearing of translated model. 
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Translated 


Extruded 


Warp Tow 



Fill Tows 



Resin 


+9.167e-01 

+8.333e-01 

+7.500e-01 

+6.667e-01 

+5.833e-01 

+5.000e-01 

+4.167e-01 

+3.333e-01 

+2.500e-01 

+1.667e-01 

+8.333e-02 



Figure 13: Comparison of predicted Tsai-Hill failure criterion of 

translated and extruded models subjected to load generated 
by 1 % £ xz shearing of translated model. 
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Chapter: 8 


MODAL TECHNIQUE FOR THREE-DIMENSIONAL 
STRESS ANALYSIS OF PLAIN WEAVE COMPOSITES 


ABSTRACT 

Textile composites present a unique challenge to the stress analyst due to their large 
microstructure. A global/local method was developed which significantly reduces the 
computational challenge and retains reasonable accuracy. This method uses homogenized 
engineering properties for the global analysis. The nodal displacements from this global analysis 
are used to determine the magnitudes of a few fundamental modes. These magnitudes are used 
to scale and superpose refined local solutions for these fundamental modes. A stubby 
cantilevered plate subjected to end moment and transverse load was analyzed using the proposed 
method. Hie results were quite encouraging. 


INTRODUCTION 

Plain weave composites consist of interlaced fiber bundles which are impregnated with a resin 
and then cured. It is not practical to discretely model each fiber bundle even when analyzing a 
simple coupon. Instead, the analysis is divided into two stages: a global analysis and a local 
analysis. The global analysis uses some form of effective constitutive properties. This might take 
the form of homogenized engineering moduli or macro elements, which are finite elements 
which approximately account for the microstructure within a single element [1,2]. To obtain 
detailed stress information requires a subsequent local analysis which discretely models the 
individual fiber bundles and matrix pockets. 

There are many global local strategies [e.g. 3-6]. A critical step in all of these strategies is how 
to apply boundary conditions on the detailed local model based on nodal forces and 
displacements from a crude global model. In an earlier study, the displacements from the global 
model were imposed along the entire global/local boundary of the local model. [3] Compatibility 
is exactly satisfied by mis procedure. However, this procedure resulted in excellent predictions 
of the stress distribirion away from the global/local boundary, but large errors near the 
boundary. This paper describes a procedure which satisfies compatibility between the global and 
local models in an average sense. In particular, the magnitudes of a few basic deformation 
modes are determined from the global analysis and then imposed on the local model. The next 
section will describe the global/local method. Then the cantilevered plate configuration selected 
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for this evaluation and the finite element models will be described. Finally, a few results will 
be presented which illustrate the advantages of the technique. 


DESCRIPTION OF GLOBAL/LOCAL METHOD 

Fig. 1 shows a schematic of the global/local procedure. The global analysis was performed using 
effective engineering properties. The effective properties for interior cells were determined using 
periodic boundary conditions for an infinite array of unit cells. [7]. Because of free-surface 
effects [3], the effective properties are different for interior and exterior cells. For the exterior 
cells the properties E,, Ey, v^, and were determined using a two unit cell model with free 
surfaces at the top and bottom (ie. z= ±2 c, where c= mat thickness). 

The nodal displacements from the global analysis were used to determine the magnitudes of the 
selected fundamental deformation modes for the region of interest. In this study the region of 
interest was simply a 1/8 subcell of a unit cell (referred to herein as simply a subcell) and was 
modeled by one of the 20-node elements in the global model. The displacements for the 20 
nodes were used to determine the 20 coefficients in a tri-quadratic polynomial fit for each 
displacement component. These polynomial curve fits were differentiated to obtain the 
magnitudes of the particular modes. In particular, for this study 12 fundamental modes were 
selected which were assumed to be constant for a unit cell. There were 6 constant strain modes: 
€ x , € y , € z , e xy , , and €„. There were also 6 -flexure modes: e xv , € xz , e yz , € yx , e zx and 

e z y . The comma in the subscript indicates differentiation. For example, the mode e xz equals 
d 2 u/dzdx evaluated at the centroid of the unit cell containing the subcell. 

Because of symmetry and the invariance in the y-direction for the configuration studied, a 
maximum of only 5 were non-zero: e x , e z , , and e^ x . The refined local model was 

subjected to unit values of each mode. For an interior cell the modal magnitudes determined 
from the global analysis were then used to scale and superpose the interior unit modal solutions 
to obtain a refined solution for the selected region. For an exterior cell this procedure was 
modified, as described later in this section. This is essentially a higher order version of the 
traditional use of micromechanics in a structural analysis. It is not unusual to obtain a global 
solution using homogenized properties and then use the calculated global strains with the actual 
material properties or a unit cell analysis to determine local stresses. The difference is that 
traditionally the unit cell analysis is performed for constant strain modes only. The advantages 
of including the linear strain modes will be discussed in the Results and Discussion section. 

It should also be noted that different boundary conditions were used to determine the interior and 
exterior unit modal solutions. The interior solutions were obtained from a model which assumed 
that the cell was buried inside an infinite array of cells. Only one exterior mode solution was 
used in this study: € x . The exterior € x mode solution was obtained using a model with two unit 
cells through the thickness (i.e. the z-direction) and an infinite number in the x- and y-directions. 
That is, there were traction-free surfaces at z= ±2 c. The possibility of more free surfaces was 
not considered in this effort. Since the exterior e x mode also contains an e z component, it was 
necessary to account for this contribution when determining the magnitude of the pure e z mode 
to superpose. If one defines to be the average Poisson’s ratio for an exterior cell, then the 
modified e z mode is (e z from the global analysis) - (e x from the global analysis). 
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CONFIGURATION 


This section describes the cantilevered plate configuration and the finite element models selected 
for this initial evaluation of the global/local method. A simple three dimensional configuration 
was selected so that reference solutions could be obtained using conventional finite element 
modeling. Fig. 2 shows the configuration considered. The cantilevered length (x-direction) was 
1.5 times the thickness. The plate was infinite in the y -direction. The nominal curvature in the 
y-direcdon was restrained to be zero by imposing zero normal displacements on the planes y= 
±A/2. The fiber tows were assumed to follow a sinusoidal path. The ratio of the tow wavelength 
to the mat thickness was 1/3. 

Two load cases were considered: constant moment and transverse end load. For both cases 
v(x, ±X/2,z) = 0. The other boundary conditions are as follows: 

Constant moment case: u(0,y,z) = 0 u(3X, y,z) = .01 z w(0,0,0) = 0 

Transverse end load: u(0,y,z) = v(0,y,z) = w(0,y,z) = 0 

Traction in the z-direction on the plane x= 3X was -16 MPa. 

Figs. 1 and 2 show the meshes used for the reference solution and the global and local models 
for the global/local analysis. All of the models used 20-node solid elements. The global/local 
analysis uses far fewer elements and degrees of freedom. The difference would be even greater 
if more refined models were used. 

The material properties used for the fiber tows and matrix pockets are shown in the following 
table. 


Tow 


Matrix Homogenized properties 



modulus 

strength" 

modulus 

strength* 

modulusfl.C.) 

modulus (E.CA 

E„ 

206.50 GPa 

1034/-689.5 MPa 

3.45 GPa 

103.4/-241.3 MPa 

29.45 GPa 

24.97 GPa 


5.17 GPa 

41.37/-117.2 MPa 

3.45 GPa 

103.4/-241.3 MPa 

29.45 GPa 

24.97 GPa 


5.17 GPa 

41.37/-117.2MPa 

3.45 GPa 

103.4/-241.3 MPa 

4.80 GPa 

4.80 GPa 


0.25 


0.35 


0.09 

0.16 

Vn 

0.25 


0.35 


0.94 

0.99 

V-a 

0.25 


0.35 


0.94 

0.99 

G« 

2.39 GPa 

68.95 MPa 

1.28 GPa 

89.6 MPa 

1.96 GPa 

1.96 GPa 

0,3 

2.39 GPa 

68.95 MPa 

1.28 GPa 

89.6 MPa 

2.64 GPa 

2.64 GPa 

G n 

2.39 GPa 

68.95 MPa 

1.28 GPa 

89.6 MPa 

2.64 GPa 

2.64 GPa 


I.C. = interior cell * Tension/compression strength 
E.C. = exterior cell 


RESULTS AND DISCUSSION 

Stress distributions were determined for a cantilevered beam using a traditional finite element 
model using the proposed global/local procedure. Two subcells were examined in detail for each 
load case. One was in the interior and one was on the exterior of the plate. These subcells are 
shown shaded in Fig. 2. The stresses were normalized by the allowable for each stress 
component and the locations and magnitudes of the largest normalized stresses were identified. 
These normalized stresses can be considered failure indices for a maximum stress failure 
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criterion. The table below summarizes the results for several of the larger values. 

Table 1. Summary of predicted critical locations. 



= actual stress component 

F.I. = failure index = actual stress/allowable stress. 

%e = error in the G/L prediction (linear inodes included) w.r.t reference solution. 

= error in the G/L prediction (only constant modes) w.r.t reference solution. 

In the table, the location is given as warp or fill tow or matrix pocket. Subsequent contour plots 
will give a more prerise indication of the location of the largest failure indices. For both load 
conditions and subcell locations the normalized <J and o 33 were largest and the largest values 
were in the warp tows. The errors were zero for the interior subcell for the constant moment 
load case. This was expected. For the transverse load case the errors tended to be significantly 
smaller for the interior subcell that for the exterior one. The most critical failure indices were 
predicted very well for the interior subcell. Although the predictions were not as good for the 
exterior subcell, the predictions were generally much better than if only constant modes had been 
used. 

Figs. 3 and 4 show normalized stress contours for the interior and exterior subcell warp tow for 
the transverse end load case. The most critical locations were in the warp tow and the critical 
stress components were <J 33 and <T 13 , so that is what is shown. The location and magnitude of 
the peak values are labeled. The figures show that the global/local procedure predicts the 
location of peak stress and the distribution quite well for the interior subcell and reasonably well 
for the exterior subcell. 


SUMMARY 

A global/local procedure was described for analyzing textile composites. The procedure would 
be equally applicable for other materials for which a characteristic unit cell can be identified. 
The initial evaluation of the method showed the technique is quite promising. It should be 
noted that in this study, a simple configuration was studied. Further evaluation is needed for 
more complex configurations. The major obstacle is obtaining reference solutions for complex 
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configurations. In fact, even for the current study only very coarse meshes were used for the 
reference and local solutions in order to limit the number of equations involved. For this reason, 
this paper did not discuss the significance of the stress components, but instead concentrated on 
showing that the proposed global/local technique predicts the same trends as a similarly refined 
traditional analysis. However, the global/local technique is far less costly than a traditional 
analysis. 
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Figure 3 Comparison of predicted normalized stresses for warp tow in 
interior subcell for transverse end load case. 


139 






140 






Chapter: 9 


Effect of Various Approximations on Predicted Progressive Failure 

in Plain Weave Composites 


Abstract 

Three-dimensional finite element analysis was used to simulate progressive failure of a plain 
weave composite subjected to in-plane extension. The loading was parallel to one of the tow 
directions. The effects of various characteristics of the finite element model on predicted 
behaviour were examined. The predicted behavior was found to be sensitive to quadrature order, 
mesh refinement and the material degradation model. Also the sensitivity of the predictions to the 
tow waviness was studied. The predicted strength decreased considerably with increased 
waviness. More numerical studies and comparisions with experimental data are needed to 
establish reliable guidelines for accurate progressive failure prediction. 

Introduction 

Textile composites consist of interlaced tows (fiber bundles) which are then impregnated with 
a matrix material and cured. Figure 1 illustrates the architecture for a plain weave composite. The 
interlacing of the tows offers the potential for increased through-thickness strength. There is also 
the potential for reduced fabrication costs, since fairly complicated shapes can be formed using 
textile machinery. One disadvantage of textiles is the difficulty in predicting their performance. 
The complex geometry makes detailed stress analysis quite challenging. The early analyses were 
based on modified laminate theory. ( e.g. References 1,2) In recent years there have been a few 
attempts to discretely model the fiber bundle architecture and predict internal stress states, (eg. 
References 3-11) Reference 1 1 presented a particularly interesting progressive failure analysis of 
a plain weave composite. The results in Reference 1 1 consisted of nominal stress-strain curves. 
The response of the composite was almost linear for in-plane extension and highly nonlinear for 
in-plane shear. The nonlinearity was primarily a result of progressive damage. However, little 
information was provided on damage evolution and load redistribution within the composite 
during the loading process. Also, there was no indication of the sensitivity of the predictions to 
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mesh refinement or other approximations inherent in such analyses. 

This paper has two objectives. The first is to evaluate the sensitivity of predicted progressive 
failure to quadn ture order, mesh refinement, and choice of material degradation model. The 
second objective is to describe the nature of the progressive failure process for two weaves with 
very different waviness. Loading consisted of a nominally uniaxial stress along one of the fiber 
tow direc tions Only mechanical loads were considered in this study. To s imp li f y the response the 
composite was assumed to consist of an infinite number of identical unit cells in all three 
coordinate directions. 

The following sections begin with a description of the basic theory used for progressive 
damage modelling. Then the configurations will be described. Finally the results from the 
numerical simulations will be discussed. 

Theory 

There is no “right” way to model damage evolution in a textile composite that is also 
practical. It is not feasible to discretely model all of the damage, so approximation is unavoidable. 
Perhaps the simplest procedure to account for damage in a finite element model is to modify the 
constitutive matrix at the quadrature points of a numerically integrated finite element. If history 
effects are not included, the analysis of the loading becomes a series of elastic analyses. Of 
course, there are many possibilities for how to modify the constitutive matrix. Three techniques 
were used herein. The first method considered the material totally failed (ie. the entire constitutive 
matrix was reduced by 3 orders of magnitude) when any allowable stress component was 
exceeded. This method will be referred to as the non-selective discount method. The second 
technique selectively reduced the rows and columns of the constitutive matrix C (where <J=C£) 
according to the particular stress allowable which was exceeded. For example, if the third stress 
component exceeded the allowable, the third row of C was set to essentially zero (reduced by 3 
orders of magnitude) to eliminate that stress component. To keep the constitutive matrix 
symmetric, the thir d column would also be set to zero. Zeroing the column has the undesirable 
side effect of stiffening the material with respect to the other stresses. This degradation method 
based on setting rows and columns of C equal to zero will be referred to as the selective RC 
method. The third technique selectively reduced the engineering moduli according to the 
particular stress allowable which was exceeded. Except as noted, this technique was used in the 
analyses. This technique will be referred to as the Blackketter method.(see Reference 1 1) 

The progressive failure analysis begins with a linear analysis of the un d a m aged configuration. 
Based on the calculated stresses, the initial load was scaled back so that failure would occur only 
at points which were within two percent of the maximum normalized stress. ( The stresses were 
normalized by the respective strengths.) The constitutive matrix was modified at the failure points. 
Residual forces were calculated and used to determine the incremental displacements. The total 
displacements were updated and used to determine the new stresses, which were then used to 
predict further failure. If no further failures occurred at the current nominal strain state, the 
nominal strain was incremented to cause failure. This procedure was repeated until the nominal 
strain exceeded one percent 

Configurations 

The fiber bundles or tows in the models were generated by translating a lenticular cross- 
section along a sinusoidal path. The waviness ratio is defined to be the ratio of the woven mat 
thickness to the wavelength. The weave consists of warp and fill tows oriented perpendicular to 
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each other. In general, the warp and fill tows could be different in terms of material and shape. In 
this study they were assumed to be identical. More details about the mesh geometry can be found 
in Reference 9. The following subsections describe the finite element meshes, the boundary 
conditions, and the material properties. 


Fini te. Element Meshes , , - . „ . t . 

Symmetry in the material and loading was exploited so that only 1/32 ofaunit cell hadtobe 

modeled. A wide range of mesh refinements were used, as shown in Figure 2. The _ciude mesh h 
only 4 elements and 42 nodes. The most refined mesh had 192 elements and 1049 nodes. The 
elements were 20-node hexahedral elements Since only 1/32 of the unit ceil was modeled, these 
refinements correspond to full cell models with 128 to 6144 elements. 


B oundary Conditions . • . 

The periodic boundary conditions for a complete unit cell are quite simple. The appropriate 

boundary conditions for a 1/32 unit cell are a bit more complicated. Derivation of the penodic 
boundary conditions is somewhat tedious, so details will not be given here Details can be found 
in Reference 9. The periodic conditions are fisted below. Figure 2 shows the coordinate system 

assumed. 


u(a/2,y,z)= Uo 
u(0,y,z) = -u(0,y,-z) 
u(x,0,z) = u(x,0,-z) 


v(x,a/2,z) = constant 
v(0,y,z) = v(0,y,-z) 

v(x,0,z) = -v(x,0,-z) 


w(x,y,c/2) = constant 
w(0,y,z) = -w(0,y,-z) 
w(0,y,z) = -w(0,y,-z) 


The load was controlled by specifying the magnitude of n. This corresponds to uniform uniaxial 
extension in the warp tow direction. 


Material Properties 

The unit cell contains two “types” of materials: the tows and the matrix pockets. Relative to 
the material coordinate system, the properties of the tows are invariant (before damage ocwr&Of 
course, the properties of the tows are needed in the global coordinate system. Fourth order tensor 
transformation formulas were used to perform the required calculations. The rotation angles to be 
used in these formulas were specified at each quadrature point by using interpolation. This 
procedure was shown in References 8 and 9 to be preferable to using a single angle for the entire 
element. The particular properties used are listed below. These properties, which are 


143 



representative of AS4/3501-6 graphite/epoxy, are from Reference 13. 


Tow properties 

Modulus Strength 


E„ 

154.27 

GPa 

2342.0 

MPa 

En 

10.80 

GPa 

56.6 

MPa 

K 

10.80 

GPa 

56.6 

MPa 

G|j 

7.47 

GPa 

48.7 

MPa 

G,j 

7.47 

GPa 

48.7 

MPa 

G* 

3.33 

GPa 

48.7 

MPa 

V|2 

0.28 




v 13 

0.28 




v 23 

0.34 





Matrix properties 
Modulus Strength 


3.45 

GPa 

84.85 

MPa 

3.45 

GPa 

84.85 

MPa 

3.45 

GPa 

84.85 

MPa 

1.28 

GPa 

101.00 

MPa 

1.28 

GPa 

101.00 

MPa 

1.28 

GPa 

101.00 

MPa 


0.35 

0.35 

0.35 


Results and Discussions 

Most of the results in this paper illustrate the effects of characteristics of the finite element 
model on the progressive failure prediction. The effects of quadrature order, mesh refinement, and 
material degradation strategy will be considered first Then the effect of the tow waviness on 
failure behavior will be discussed. Except where indicated otherwise, the results will be presented 
for a waviness ratio of 1/3 and for the material degradation strategy in Reference 1 1, ( the 
“Blackketter method” ). 

Figure 3 shows the effect of quadrature order on the stress-strain curve. The peak stress 
obtained using 8 quadrature points (2x2x2) is 3 percent higher than that obtained using 64 points. 
The peak stresses obtained using 27 and 64 quadrature points differ by 1 percent Damage 
initiati on is predicted 3 percent earlier when 64 points are used. After the large stiffness loss 
which occurs at about 0.3 percent strain, there are even larger differences in the predictions. In 
Reference 12, non-selective discount was used for the same configuration. The difference in the 
peak stress obtained using 8 and 64 point quadratures was 10 percent Hence, the sensitivity of the 
predictions to quadrature depends on the degradation model. The sensitivity to quadrature order 
is not suprising. For example, when more quadrature points are used, the more extensive 
samp lin g is more likely to find the extremes in the stress field. One might expect a refined mesh 
to exhibit less sensitivity to quadrature order than a coarse mesh. For the meshes considered in 
this study, this was the case. Also, when failure occurs within an element, and the constitutive 
matrix is modified, the element becomes inhomogeneous. The numerical integration effectively 
fits a polynomial function to the variation of the material properties. Since the material properties 
are very different in the failed and unfailed parts of the element, it is difficult to obtain a good fit 
In fact, there is a concern as to whether the assumed quadratic displacement functions for a 20- 
node element are sufficient to obtain a reasonable approximation regardless of the integration 
order. 

Figure 4 shows the effect of mesh refinement on the predicted stress-strain curve for the two 
waviness ratios. The four-element model predicts the correct trends, out is quite inaccurate. The 
error in the peak stress is much worse for larger waviness ratio. The peak stress and 
corresponding strain for a coarse mesh is larger than for a refined mesh, but the peak stress and 
corresponding strain for a moderately refined mesh are not necessarily bound by the extreme 
cases of mesh refinement, (see Figure 4b) The peak stress and the corresponding strain decrease 
with increased waviness. The ratio of the initial damage stress to the peak stress and the 
corresponding ratio of strains increase with decrease in waviness ratio. That is. not only is the 
strength reduced by increased waviness, but there is also earlier damage initiation relative to the 
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peak stress. After the peak stress is reached, there is a precipitious drop in stress with little change 
in st rain For the 1/3 waviness ratio and a coarse mesh, this “collapse” actually occured in stages, 
with a s mall strain increase between the stages. For the most refined model the collapse occurred 
in a single stage. For the 1/6 waviness ratio the collapse occured in a single stage for all three 
mesh refinements. It should be noted that in real specimens there would be a distribution of 
strength in the various cells. Accordingly, the nominal stress-strain curves would not be expected 
to have such deep jogs. Also, it would not be possible to impose displacement or load control for 
the individual cells. Load redistribution and the stability of the damage growth would determine 
whether the nominal stress-strain curve would appear brittle or there would be significant 
nonlinearity before total fracture. 

Figure 5 shows damage volume versus the nominal strain for the warp and fill tows for four 
mesh refinements. Damage volume for the resin is not shown since there is almost no damage. 
The curves are quite close for the 108 and 192 element meshes, which suggests that the curves 
might be close to convergence. As the mesh is refined, the increments in damage volume become 
smaller , but more numerous. The damage volume just after collapse tends to decrease with mesh 
refinement for both warp and fill tows. This is not surprising since the corresponding strain is also 
smaller. The curve for the 192 element mesh in Figure 5a has points labeled A and B. From A to 
B the damage initiation and growth was dominated by the inter-tow normal stress O33. The sudden 
increase in damage at point B was due to the stress component a ( 3 ( the x 1 direction is parallel to 
the fibers). 

The damage volumes just after the peak stress is reached is as follows 

Waving Ratio Warp Fill Percent Stiffness Loss 

1/6 -76 .88 41 

1/3 .41 .21 45 

Interestingly, the damage volume is larger for the 1/6 waviness ratio, but the percentage stiffness 
loss is less. 

Figures 6 and 7 show the effect of mesh refinement and waviness ratio on damage 
accumulation during loading. To plot the damage zones accurately, each element is sub-divided 
into 27 blocks, where each block represents the volume associated with a quadrature point. The 
black region indicates the damage zone. The stress-strain curve for a particular mesh is shown 
above that mesh. The points labeled A,B and C indicate the correspondence between the strain 
level and the damage zone. The damage zone corresponding to point A indicates the initial 
damage. The damage zones corresponding to points B and C indicate the pre- and post- collapse 
damage states. Also indicated are the stress components which contributed to the damage. The 
four-element mesh does not model the initial failure well for the two waviness ratios of 1/3 and 1/ 
6. The 32 element model performs reasonably well for obtaining qualitative results. Further 
numerical studies are needed to determine how close the 192 element results are to convergence. 
For the 1/3 waviness ratio C33 dominates the initial failure. This initial failure appears to be an 
inter-tow failure resembling delamination in laminated composites. The collapse is characterized 
by Cj3 failure in the warp tows. For the 1/6 waviness ratio, die collapse is due to a significant 
failure of the warp tow due to and cracking of the fill tow due to 022- 

Figure 8 shows the stress-strain curves obtained using the three degradation models described 
earlier non-selective, selective RC, and the Blackketter method. The non-selective method 
predicts 22% lower peak stress than the Blackketter method. Selective reduction of rows and 


145 



columns in the constitutive matrix results in a much larger residual stiffness just after the peak 
stress than the other two methods. Selective reduction of rows and columns in the constitutive 
matrix also does not result in a large sudden drop in the stress after the peak stress is reached. 
Concomitantly, for the selective RC method there is also no sodden increase in the damage 
volume after the peak stress is reached. (This is not shown on the plot) 

Concluding Remarks 

Simulation of progressive failure in a plain weave composite is extremely complex. 
Consequently, only approximate treatment is practical at this time. One of the goals of this paper 
was to ex amin e the effect of several approximations on predicted behaviour. One obvious 
conclusion from this study is that the predictions are quite sensitive to a number of decisions 
which must be made when assembling a finite element model. Further numerical experiments and 
comparisons with experimental data are needed to establish guidelines for accurate analysis of 
progressive failure. 

Another objective of this paper was to describe the effect of tow waviness on damage 
accumulation. The results suggest that the degree of waviness not only affects the stress at which 
damage initiates, but also the type of damage which occurs. Also, the stress component 
responsible for damage changed during the progressive failure process. 
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Figure 3 Effect of quadrature order on nominal stress-strain curve. 

Waviness ratio ■ 1/3. Blackketter discount method was used. 
The model had 192 elements and 1049 nodes. 
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Number of elements = 4 Number of elements = 32 Number of elements = 192 

Figure 7 Effect of mesh refinement on damage development for waviness ratio = 1/6 
Blnckkcltcr discount method was used. 
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Part II 




Chapter: 1 

Part II Software Documentation 


Stress Analysis of Woven Composites 

(SAWC) 


Introduction and Abstract 

SAWC refers to a collection of three programs which were developed to the 

. f th failure behaviour of plain weave composites. The programs are a mesh gener 
study of the failure oenaviour u y r eranhical pre- and postprocessor 

(PWMeshGen), a finite element program (Hex94), and a grapmc<u pic y y 

(P1 °' 95> The primary capability of the software is the generation of solid models of plain weave 
from afew input parameters, the analysis of the solid model, and the display of finite 
element meshes before and after deformation and stress distributions. Some first order progres- 
s Sure capability is also implemented. The finite element code uses a displacement fotmula- 
Z TSes rioparametric 20-node elements and social "macro” elements which can 

aCCOUn, The “r rresV;t“ form, which is general antomaticaUy. The 
user simply fiUs in a few blanks in the fonn with any text editor and then runs the program. The 
finite element program requires a complex collection of input data. Fortunately, the mesh genera 
tor* provides most of the data. The plotting program was originaUy devetoped tismg ^D paPhtcs 
Ubrarv and a painter’s algorithm to manage hidden fine removal. Hardcopy of the displ y 
obtained by writing a PostScript file, which is then printed. Although the program now uses the 
rT cmanhics library (which has built-in hidden line removal), 2D primitives are still used. 

» a combination of C. C++ , and Foriran 77. The pro^ 
were develoLonan IBM RS/6000 using the ADC (IBM UNIX) operating system version 3.2.5. 
However PWMeshGen and Flex94 should port easily to any Unix workstation The plotting pr 
gram, Plot95, uses the graphics library openGL and the Motif user interface. The openGL library 

and Motif user interface is available on major workstations. 

This software can be installed on a workststion without root (super user) access. 

This manual contains four sections: one section for each program and then a final section 
which describes installation of the software and sample files. 

Listed below are papers which describe the theory and example application of the software. 
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Notice of Restrictions 

The Motif interface was developed using AlXwindows Interface Composer. This pro- 
gram generates Motif code which is than compiled and linked with the other parts of the Plot95 
program. Appropriate copyright notices for the software generated are shown below. 

/* 

* COMPONENT_NAME: AIC AlXwindows Interface Composer 

* 

* ORIGINS: 58 

* 

* 

* Copyright IBM Corporation 1991, 1993 

* 

* All Rights Reserved 

* 

* Permission to use, copy, modify, and distribute this software and its 

* documentation for any purpose and without fee is hereby granted, 

* provided that the above copyright notice appear in all copies and that 

* both that copyright notice and this permission notice appear in 

* supporting documentation, and that the name of IBM not be 
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* used in advertising or publicity pertaining to distribution of the 

* software without specific, written prior permission. 

* 


IBM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING 
ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 
PURPOSE. IN NO EVENT SHALL IBM BE LIABLE FOR ANY SPECIAL, INDIRECT OR 
CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF 
USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR 
OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE 
OR PERFORMANCE OF THIS SOFTWARE. 


*/ 

/* 

* $Date: 93/04/26 21:05:25 $ SRevision: 2.19 $ 

* 

* 

* Copyright (c) 1992, Visual Edge Software Ltd. 

* 

* ALL RIGHTS RESERVED. Permission to use, copy, modify, and 

* distribute this software and its documentation for any purpose 

* and without fee is hereby granted, provided that the above 
♦copyright notice appear in all copies and that both that 

* copyright notice and this permission notice appear in supporting 

* documentation, and that the name of Visual Edge Software not be 

* used in advertising or publicity pertaining to distribution of 

* the software without specific, written prior permission. The year 

* included in the notice is the year of the creation of the work. 



Chapter: 2 

User’s Manual for PWMeshGen 
(Plain Weave Mesh Generator) 

A collection of Fortran and C programs were developed to expedite the generation of finite 
element meshes for plain weave composites. These programs are currently intended to be run 
under the UNIX operating system. However, only a few changes are required for other operating 
systems. 

The tow path of the plain weave is assumed to be sinusoidal. The user can select between 
translated and extruded tows (see Ref. 1). Although several programs are required to generate a 
mesh, an executive program has been provided which orchestrates the transfer of data from one 
program to the next Hence, the collection of programs appears to the user to be a single program. 
The executive is named PWMeshGen. 

To simplify program operation, the input file is a form. This form contains labels which 
remind the user of the required data and order of the data. To obtain a copy of the form, simply 
execute the program "PWFonn” with the command line parameter filename , which will be the 
name of the generated form. For example, executing the command 

PWFonn inFile 

will generate a file named <inFile>, which is shown in Figure 1. Figure 2 shows a typical finite 
ftle.nne.nt mesh with labels which define the terms in the generated form. 

When this file is edited, it is critical than none of the form labels be changed, since the 
labels are used to guide input Once the form is complete the mesh is generated by executing the 
command 

PWMeshGen inFile 

The programs will generate several files, which will be discussed in the next section. Sample data 
in the sub-directories “Sample3” and “Sample4” on the distribution media include completed 
forms. The file names for these forms are Samples/Sample3/Input/meshl and S amples/S ample4/ 
Input/mesh2. 

Warning 

Early in the mesh generation process there are duplicate node numbers. One of the tools 
removes the duplicate node numbers. This tool uses a tolerance to determine whether two points 
are coincident. This tolerance is hardwired to be .00001. It can be changed by editing the file 
MeshClass.C in PWMeshGen/MeshToolSource. Line 140 is 
#define EPSILON le-5 

To change the tolerance, simply change this value and recompile. 

Output Files: 

Several output files are generated. These files are for use with the finite element program 
Flex94. The following files are generated during a typical execution: 
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File 

new.flex 
new.sflx 
new. as 
new. am 
mat_list 
new.flx 


xExtension.mpc 

xyShear.mpc 

xzShear.mpc 

ExtSingleConstraints 

xyShearConstraints 

xzShearConstraints 

xExtension .Loads 

xyShear.Loads 

xzShear.Loads 

eighth.flx 

eighth. am 

eighth.as 


Description 

Main input file for Flex94. 

Mesh file for Flex94. 

Element rotation angle file: single angle. 

Element rotation angle file: multiple angle. 

A material list of the elements. 

A simple mesh file used for plotting the mesh 
and determining boundary conditions automatically. 

Multipoint constraints for extension in the x direction. 
Multipoint constraints for in-plane shearing. 
Multipoint constraints for transverse shearing. 

Constraints for extension. 

Constraints for in-plane shearing. 

Constraints for transverse shearing. 

Loads for extension in Ac x-direction. 

Loads for in-plane shearing. 

Loads for transverse shearing. 

A simple mesh file for plotting the l/8th unit cell. 

Element rotation angle file for l/8th unit cell: 
multiple angle. 

Element rotation angle file for l/8th unit cell: 
single angle. 


References: 

1 Chapman C 1993. Effects of assumed tow architecture on the predicted moduli and 

stresses in woven composites. Master thesis, Departrae* of Aerospace Engineering, Texas 
A&M University. 
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Input File for Mesh Generation Program 


Thickness of mat: 

Waviness ratio: 

Tow type: 

Tow elements in z- direct ion: 

Primary elements in y-direction: 

Resin elements above and below tows: 

Execution flow flags: Type yes beside functions to be performed. 
Generate 1/8 unit cell: 

Renumber nodes to reduce profile of stiffness matrix: 


Notes 

Tow type: 1=> extruded 

2=> translated 


★*****★*♦★** 


* *End of Input File For Mesh Generator 




Figure 1 : Form used to define input for mesh generator. 
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Appendix A: Use of output files with Flex94 


The files generated using PWMeshGen are used in conjunction with the finite element 
program Flex94 and the mesh plotting program Plot95. The file new.flex is the main input file for 
Flex94. An example of this file is shown in Figure Al. This example specifies that the mesh 
new.sflx will be subjected to extension in the x-direction as indicated by lines 6, 26, and 29. Mod- 
ifications for other load cases are given below: 

For in-plane shear, these lines would need to be changed as follows: 

6 ‘xyShearConstraints’ 

26 ‘xyShear.mpc’ 

29 ‘xyShear.Loads’ 

For transverse shearing, these lines would need to be changed to: 

6 ‘xzShearConstraints’ 

26 ‘xzShear.mpc’ 

29 ‘xzShear.Loads’ 

And finally, for extension in the z direction, the lines would be: 

6 ‘ExtSingleConstraints’ 

26 ‘zExtension.mpc’ 

29 ‘zExtension.Loads’ 

Note that ExtSingleConstraints is also used for extension in the x-direction. 


165 



1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 


4 3 

'Title Here' 

' altemate_input * 

' new. sf lx ' 

' altemate.input ' 

• ExtSingleConstraints ' 

' end_mesh_ input ' 

• 3d ' 

1 

206. 9e9 5.171e9 5.171e9 .25 .25 .25 2.386e9 2.386e9 2.386e9 0000000 
2 

5.171e9 206. 9e9 5.171d9 .00625 .25 .25 2.386e9 2.386e9 2.386e9 0 0 0 0 0 

0 0 
3 

3 45e9 3 . 45e9 3.45e9 .35 .35 .35 1.28e9 1.28e9 1.28e9 0000000 
0 ' end:T300/5208 1 

1 end_material_input ' 

' loop ' 

1111 
2 2 2 1 
3 3 4 1 
0 0 0 0 
' end_pick* 

' altemate_input 1 
’ new . am ' 

• altemate_input 1 
' xExtension . mpc ' 

• end_of_misc_options ' 

1 alternate_input ' 

1 xExtension . Loads * 

' end_loads ' 

• end ' 


Figure Al: Typical new.flex file generated with PWMeshGen. 


166 



Chapter. 3 

FLEX94 

User’s Manual 


Command: fe size 

Comments: size = maximum number of terms in global stiffness matrix. If size is omitted, a 
default size is assigned by the program. The default size is 1500000. 

The analysis of an infinite array of unit cells only requires a single mesh. Such analysis is 
useful for determining homogenized engineering properties and stress (or strain concentrations). 
Such analysis proceeds much like traditional finite element analyses except that the boundary 
conditions are fairly complicated. Utilities have been developed to automatically generate the re- 
quired boundary conditions for various load conditions. 

Global/local techniques were developed as part of this research project. There are many 
possible global/local methods. The ones evaluated used macro elements (Refs. 1-3) in the global 
mesh and ordinary finite elements in the local mesh. Two types of macro elements are supported: 
single field ( Refs 1, 2) and multi-field (Refs 3). After a global analysis is performed using macro 
elements, the detailed stress distributions within a weave unit cell are determined using displace- 
ments or forces from the global analysis to determine the boundary conditions for the local mod- 
els, which include details of the weave architecture. The global/local analysis software was not 
sufficiently automated to release it as part of this software. However, the macro elements are in- 
cluded in the finite element program. Reference 4 discusses one of the more promising proce- 
dures evaluated. 

The input file for Flex94 can be broken into several blocks which must appear in the fol- 
lowing order: 

1. Mesh Input 

2. Material Properties 

3. Macro Element Input (Optional) 

4. Miscellaneous Options 

5. Loads 

6. Macro Element Data 

7. Failure Analysis 


A description of each block is given in the following sections with an example. 
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1. Mesh Input 


Example: 

‘analysis_type’ 

‘LINEAR’ 

‘end_options’ 

43 2 

’Small Mesh’ 

’attematejnput’ 

’sml.msh’ 


Example; 

’standardjnput’ 


’coordinates’ 
102 2 

1 -1.5 -1.5 

2 -1.5 -1.3 

3 -1.3 -1.5 

4 -1.1 -1.5 

101 1.5 1.3 

102 1.5 1.5 


’connectivity’ 

1 4 1234 

2 4 2336 
344876 

42 4 97 95 99 100 

43 4 100 99 99 101 


’define_element_type’ 
243 1 43 1 

0 0 0 0 


Description: 

This option allows the user to define the analysis type. 

‘LINEAR’ selects linear analysis. This has to be replaced 
by ‘SELECT for Selective discount method or by ‘NSELECT 
for Non-Selective discount method. 

NumberOfElements DegreesOfFreedomPerNode 
Title - Must be in single quotes. 

This option allows user to put mesh in another file. Filename must 
be on the following line. At end of file ’sml.msh’, ’standardjnput’ 
should be used to return to original file, ’alternate Jnpir: can only 
be used in the original Flex94 input file (eg. you could not use the 
command in ’sml.msh’) 

Description: 

If ’alternatejnpuf was used, ’standardjnput’ would return input to the 
original input file, ’standard Jnpuf can only occur where a command is 
appropriate. For example, it could not appear in the middle of reading 
coordinates. 

Command to signal start of coordinates. 

NumberOfNodes NumberOfCoordinateDimensions 
NodeNumber Coordinates 


Command to signal start of connectivity. 

Element# NumberOfNodesInElement Connectivity 
Connectivity must be specified in clockwise order for 2D elements. For 
20-node 3D elements, the order of the nodes is shown in Figure 1 . 


Command to start element type definition. 

ElementType FirstElement LastElement Increment (In this case, ele- 
ments 1 throught 43 are of type 243). 

End with four zeroes. The relevant element types are Isted below: 

243 2D element 

300 3D element 

851-899 single field 

801-849 multi-field 
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’select_quadrature_order’ Command to start selection of quadrature order for each element. 
2 i 43 i QuadratureOrder FirstElement LastElement Increment 

0 0 0 0 End with four zeros. To obtain stresses, a quadrat jre order of 

2 for 2D and 3 for 3D has to be used. 


’single_const paints’ 

40 1 0 

41 1 1 
18 0 1 
102 1 1 
100 0 1 
101 1 0 
0 0 0 


Command to set single constraints on individual nodes. 
NodeToConstrain ConstraintDofl ConstraintDof2 ... 

1 == Constrain Dof 
0 = Don’t Constrain Dof 

Note: The number of constraints at each node must be equal to the 
number of Dof per node which was set at the beginning of the 
mesh block. Example shown is for 2 Dof per node. 

End with zeros. 


’plane’ 

1 -1.5 1 

1 -1.5 2 

2 -1.5 2 
2 -1.5 1 
0 0 0 


Command to set constraints on a plane, 
idir coord jeon 

idir = direction of normal to plane in (xl ,x2,x3) space 
coord = coordinate of plane 
jeon = restraint direction 
end with zeros 


’end_meshjnput’ exit this input section 


2. Material properties 

This section defines the material library and which elements have which material proper- 
ties. Flex94 was designed to handle various types of constitutive definitions (eg. 2D, 3D, proper- 
ties for a beam, etc.) However, for textile analysis only one option is relevant - ’3D’. This option 
requires the 3D elastic properties to be given as shown below. For 2D analysis the 3D properties 
which are input are used to determine the 2D properties for plane strain analysis. 


Example; 


Description: 


’3D’ 

i 

206.9e9 5.171 e9 5.171e9 
.25 .25 .25 

2.386e9 2.386e9 2.386e9 
0 

000 

000 


Command to start reading of 3D material properties. 

Material group number used later in assigning properties to elements. 
Young’s Moduli (E n E22 E33) 

Poisson’s Ratios ( v 12 v 13 v 23 ) 

Shear Moduli (G 12 ®i 3 G 23 ) 

Rotation about z-axis (z-axis is out of plane for 2D problems) 

(thermal expansion coefficients.. .not used or implemented) 

(moisture expansion coefficients...not used or implemented) 


2 Next material group 
5.171e9 206.9e9 5.171e9 

.00625 .25 .25 

2.386e9 2.386e9 2.386e9 

0000000 

3 

3.45e9 3.45e9 3.45e9 
.35 .35 .35 


169 



1 .28e9 1 .28e9 1 .28e9 
0000000 
0 

’end_materiaLinput’ 

’loop’ 

31 431 
1 2 43 4 
1 3 43 4 
0 0 0 0 
’end_pick’ 


Give zero as material group number to end input. 

End input of material properties 
Command to start specifing material group. 
MaterialGroupNumber FirstElement LastElement Increment 


End with zeroes 

End selection of material properties for elements 


Comments: For a mesh consisting of macro elements only, there is no need to input material 
properties. (It will do no harm, but the data will not be used.) Hence, the following lines are suf- 
ficient for the material property section. 


’end_materialjnput’ 
’end _pick* 


3. Macro Element Input 

Most of the data for macro elements will be specified in another file, as described shortly. 
The following must be included in the main input file if macro elements are being used. 


Description: 

Command to start reading of macro element mesh, 
macro element type 851-899: single field 
801-849: multi-field 

NumberOfNodes NumberOfElements NumberOfDimensions 
NumberOfDofPerNode 

length of connectivity array length of coordinate array 
Minimum requirements are: 

Connectivity: numberOfElements * (numberOfNodesPerElement + 9)+1 
Coordinates: numberOfNodes * numberOfDimensions + 2 

number of elements in macro element submesh 
number of degrees of freedom per node in macro element submesh 

’title' 

’alternatejnput’ 

’name’ name of alternate input file 

(what is in this file will be described in section “6”) 

Repeat above commands of section 3 for each type of macro element 
to be used. 


F.x ample: 

’read_macro_mesh’ 

851 

103 95 2 
2 

2500 500 
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initmacro’ 

2 

1 2 


NumberOfMacroElementTypes 

List of elements which need to be initialized 


4. Miscellaneous Options 


Element Material Rotation Angle: For the analysis of textile composites, the material proper- 
ties of the elements making up the tow are the same in the material coordinate system. These 
properties must be transformed to the global coordinate system. Flex94 allows the user to specify 
the angular orientation of the elements. For 2D, the user can specify the angle of rotation for an 
entire element only. For 3D, however, Flex94 also allows the user to specify the angle of rotation 
for each node in an element. The angle of rotation may be specified using three different com- 
mands: ’angles2d\ ’angles3d\ and ’angles_multiple\ The angles are specified in terms of de- 
grees. 


Example; 

angles2d’ 

1 0.00 

2 5.17892 

3 10.28684 

4 15.34983 


Description; 

Command allows the user to specify the angles for a 2d analysis. When 
using this option, angles specify the rotation about the z-axis. 
(Out of plane.) Angles must be specified for all elements in the 
mesh and are positive for a dock-wise rotation. 

ElementNumber Rotation AboutZAxis 


42 -5.17892 

43 0.000 


Example; 

’angles3d’ 

1 1 0.00 

2 2 5.857 

3 3 6.449 


Description: 

Command allows the user to specify the angle and axis of rotation for 3D 
analysis. Again, the angles must be specified for all elements. 
ElementNumber AxisOfRotation Angle 


42 1 0.00 

43 2 2.48 

'angles_multiple’ 
1 2 20 
6.724670 
7.294361 

4.009413 
0.000000 
2 1 20 


Command allows the user to specify the angles of rotation for 3d. 
ElementNumber AxisOfRotation NumberOfAnglesForElement 
Angle(l) 

Angle(2) 

Angle(19) 

Angle(20) 

Angle(n) corresponds to the rotation at the nth node specified 
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5.877652 

2.332992 


in the connectivity of the element. 


42 1 20 0.000 0.000 0 .... 

43 2 20 2.489 2.476 2.4... 

It is often more convienient, when specifing the material rotation angles for elements, to 
use ’altemalte_input’ to allow the angles to be kept in another file. When doing this, remember to 
put ’standardjnput’ at the end of the file to let Flex94 return to the original input file. 

Multipoint Constraints: Another miscellaneous option which Flex94 allows, is the specification 
of multipoint constraints. When specifying multipoint constraints, the user must specify a master 
node, slave node, the particular degree of freedom (dof) to constrain, and a difference between the 
two dof ’s. The particular dof being constrained (ie. the slave node) cannot have been previously 

constrained. 


It is also possible to apply a mpc such that the displacement of the slave node dof is the 
opposite that of the master node dof. This is done by putting a minus sign in front of the master 
node as shown in the following example. 

Description: 

Command to start reading of multipoint constraints. 

SlaveNode MastertMode DofToConstrain Difference 

This line constrains Node 4 dof 2 to the negative displacement of Node 
1 dof 2 plus a difference of 0.1 50 

Use four zeros to signal end of multipoint constraints. 


Example: 

’mpc’ 

2 1 1 0.000 
3 1 1 0.000 
4-12 0.150 
100 -1 1 0.000 
101 100 1 0.000 
102 -1 2 0.150 
0000 


Ending Miscellaneous Options: This command must appear at the end of the Miscellaneaous 
Options section. It is shown below. 

Example: Description: 

’end_of_misc_options' Command to end Miscellaneous Options. (NOT OPTIONAL!) 

As stated earlier, it may be more convienient to keep sections of miscellaneous options in 
another file. This can be done using ’altematejnput’ with ‘standardjnput’ as explained in sec- 
tion 1. 


5 . Loads 

Various types of loads can be applied with Flex94. Some of these include the specification 
of nodal displacements and point forces. All the command options in this section are optional. 
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’altematejnput’ may be used at any time where a command can be accepted. Remember to re- 
turn to the original input file with ’standarc_input\ 


Point Forces: Point forces allow the user to specify the nodal force at a node. 


Example: 

’poinf 
1 le7 1 
3 2.345e6 2 

87 6.456e8 1 
0 0 0 


Description: 

Command to start reading of point forces. 
NodeNumbe* Force DofNumberForNode 


End reading rf point forces with three zeros. 


Displacements: Displacements may also re specified at specific degrees of freedom. In the input 
of the mesh in section 1, constraints can be input. This reduces the actual size of the problem. 
Specified non-zero displacements are also a type of constraint, but in order to reduce the problem 
size, the dof must be constrained in the mesh section also. 


Example: 

’displacement’ 
1 3.13e-3 2 

1 .025 1 

2 0.56e-2 2 
87 0.13e-2 1 
102 0.13e-2 1 
000 


Description: 

Command tc start reading of displacements. 
NodeNumbe* Displacement DofNumberForNode 


End reading rf displacements with three zeros. 


Plane Displacements: Displacements mzy be applied to an entire plane in a particular direction. 
This is known as a plane displacement This option works in conjunction with setting plane con- 
straints in the mesh input section. 


Example: 

’planeDisplacement’ 
1 -1.5 .015 2 

1 -1.5 .010 1 

2 -2.0 -.013 2 
0000 


Description: 

Command tc start reading of base displacements. 

CoordinateN-mber CoordinateValue Displacement Direction 
<_This line rdicates that on the plane x=-1 .5 specify a displacement of 
0.01. in the x-direction. 

End reading rf base displacements with four zeros. 


Linearly Varying Displacements on a Plane: Displacements may be applied to an entire plane 
so that the variation of the specified displacements changes linearly with the value of the coordi- 
nates which are parallel to the plane. For example, one may want to specify an x displacement on 
a plane x=1.5 which varies linearly with y. Displacements are calculated as d, = a y, + b where a 
and b are specified by the user and d* and y, are the calculated displacement and y coordinate at a 
specific node on the x=l .5 plane. 


Example: 


Description: 



MinearPlaneDisplacement' 

1 1.5 1 2 .1 - 05 

2 1.5 1 1 01 -.01 

000000 


Command to start reading of linearly varying plane displacements. 
This line specifies that on the plane x^l.5, a displacement in the x, 
direction given by d, = .1x 2i - .05 is being specified at each 

node i on the plane. . . 

End reading of linearly varying plane displacements with six zeros. 


To end reading of loads, ’endjoads’ must be at the end of the loads section. 


6. Seperate Input File For Macro Element Data 

Much of this file is identical to the sections described above. Hence, references will be made to 
the sections above rather than repeating all of the details. 

Mesh input block refer to Section 1: 

Comments: 

1. Do not input any restraint information. 

2. The nodal coordinates must be noimaized coordinates (eg. they must range 
between -1 and 1.) 


Material properties block.......refer to Section 2: 

numberOfNodesInMacroElement: The number of nodes in the macro element must be 
specified. It is not the number of nodes in the submesh. 

Miscellaneous options block.......refer to Section 4: 

Comments: 

1. The material rotation angles for rhe elements in the submesh is input in this 
section. 

2. Do not apply multipoint constraints to a macro element mesh. 


7. Failure Analysis 

This section describes the data required for progressive failure analysis. 

Ac ^ 1 ^ in ‘Mesh Input section 1,’ the analysis type ‘LINEAR’ has to be replaced with 

^ her ‘SELECT’ or ‘NSELECT’ option. The option ‘SELECT’ represents die ‘selective discoun 
... d ‘NSELECT’ represents the ‘Non-Selective discount method. One additional input 
fi,:“ q :luln^J‘s«re„gfhd a ta>. conrams a Us, of sirengdr values for each of rhe 

material groups used. 


174 



Example: 


Description: 


3 NumberOfMaterialGroups 

500 50 50 60 60 60 (tensile strength) aii,a 22 »<* 33 *( shear strength) <Ji 2 ,O 13 , 023 , 

-500 -50 -50 (compressive strength)a y 1 , 022.^33 


When progressive failure analysis is performed, the following additional files are created, 
‘stressstrain’ : Data file used to plot ‘nominal stress vs nominal strain’ curve. 


Example: 


Description: 


1 0.0e6 0.00 ReferenceNumber, NominalStress Value, 

2 1.3e6 0.10 NominalStrainValue (percent) 


‘ damagefield ’ : Damage progression sequence is recorded. This file may be used to study the fail- 
ure mechanism and used for graphical simulation of failure progression. 


Example: Description: 


1 1 

001 043 

000000 

650000 


ElementNumber, MateriaIG roupNumber 

Each row represents an integration point of the element. Each column 
represents a stress component. <J 11 ,o 2 2 .a 33 .Oi 2 ' a i 3 > 0r 23 is the order 
the stress components for each row. The numbers 1,4, and 3 
correspond to the first , fourth and third points on the stress-strain curve. 


2 1 

870000 

‘fcontour.n’ n = 1,2,..., number of points on the stress-strain curve : This file 
contains the contour data required to plot failure contours for each point on the stress-strain curve. 
The file format is the same as the stress contours file ‘stress’. 
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Chapter: 4 

Plot95 


User’s Manual 


Executable: plot95 

The program Mesh.app was used in debugging finite element meshes and postprocessing the 
results from finite element analyses. The program was developed for the IBM RS/6000 workst- 
stion using the Motif interface and the openGL graphics library. The program should compile 
and function properly on other workstations which use Motif and which have the openGL graph- 
ics library. However, this is would need to be verified for any particular workstation. 

The primary functions of the plotting program are: 

1. Plot a finite element mesh. 

2. Plot a deformed finite element mesh (ie. with scaled nodal displacements added to 
original nodal coordinates.) 

3. Plot stress contour lines. 

4. Plot stress contour bands. 

This program was designed to work with the finite element program “Flex94”. In brief, the 
current version of the plotting program supports the following: 

Mesh plotting for the following elements 
truss 
frame 

triangular and quadrilateral 2D elements with any number of nodes 
20-node hexahedral elements 

Contour plotting for the following elements 
4 and 8 node quadrilateral elements 
20-node hexagonal elements 

The following pages describe the use of Plot95. There are two aspects to using the program: 
preparation of the input files and interaction with the graphical user interface to obtain the type of 
plot desired. This manual will begin with a discussion of the input files followed by a description 
of the graphical user interface (GUI). 
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Input Files 


This section describes the following types of input files: 

1. Mesh file 

2. Nodal displacement file 

3. Contour data file 

In each case a fragment of a typical input file will be listed and explained. The actual input data is 
in small type and the comments are in italics. 


Mesh file: 

1941 384 3 

1 -.75000 -.65625 -.5000 

2 -.65625 -.75000 -.5000 

3 -.75000 -.75000 -.5000 

4 -.75000 -.65625 -.4951 

5 -.75000 -.75000 -.3750 

6 -.65625 -.75000 -.4951 


numberOfNodes numberOfElements numberOfDimensions 
node # coordinates (for 2D only xy coordinates are needed. ) 


element # numberOfNodesPerElemertt connectivity 

1 20 3 4 11 13 35 32 31 5 6 14 38 34 23 24 27 29 43 40 39 25 

2 2031 32 35 37 112 110 108 33 34 38 111 107 39 40 43 45 120 118 116 41 

3 20 II 12 5 53 63 36 35 13 14 54 66 38 27 28 59 61 66 44 43 29 

4 2035 36 63 65 124 114 112 37 38 66 123 111 43 44 67 69 128 122 12 45 

5 2051 52 147 149 159 64 63 53 54 150 ... 


Optional input for mesh file- 

world 
0.000.00 
1.500 1.500 


inactive 
1 64 1 
129 3S4 1 
000 


input range of screen coordinate system 
lowerLeftX lower LeftY 
upperRighiX upper RightY 

(if left out , program will automatically pick coordinates) 

deactivate elements ( elements will not be plotted) 
first last increment 
next group to deactivate 
end option with three zeros 


active 
24 321 
000 

set colors 
5 1 32 1 


0 0 0 0 

elementNodalValues 

set. values 
3 2 

1 2.334 4.566 1.13e9 

2 2.567 4.877 I .15e9 


activate elements 


Set Element Colors 
co lor Index first last increment 
(0<= colorlndex<l2) 

colorlndex corresponds to material group number 
end with four zeros 

input contour data (see format in Contour Data File Section) 

Sets values for each element Use Shade Elements to view. 
numberOfColumnsOfData selectedColumn 
total # of columns = numberOfColumnsOfData + 1 

If selectedColumn<0 t absolute value of data is input . 
Range is selected automatically. 
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ft*. values 
3 2 

0 10.5 

1 2.334 4.566 1.1 3e9 

2 2.567 4.877 I.l5e9 


Fix Values for each element . Use Shade Elements to view. 
numberOfColumnsOfData selectedColumn 
min Value maxValue 

If selectedColumn<0, absolute value of data is input 


Nodal displacement file: 

1 -.79272E-24 .13041E-02 - 42SI0E-22 nodeNum (u,v,[wj) displacement 

2 87106E-02 -.97213E-24 -.43845E-22 (w displacement is optional in 2D) 

3 -.I2200E-22 -.24332E-23 .10482E-22 

4 .37197E-22 .1 1482E-02 -.30051E-02 

5 .42536E-22 13039E-22 -.23021E431 

6 83812E-02 11663E-22 -.36294E-02 

7 .19218E-24 .5432IE-02 13092E-22 

8 15420E-23 .10875E-0I -.12807E-22 
9-.14467E-23 .45672E-02 - 52301E-02 
10 .11192E-01 63668E-02 -.49632E-22 

Contour data file: 

This section may be included in the mesh file or as a stand alone file. To include this in the mesh 
file, the option elementNodalValues must be used. 

(etc meotNodal Values) Only include if in the mesh file. 

3 Number of columns 

l Column to be input 

fixed (These 2 lines explained in Scaling options.) 

-4e7 4c7 

l j elementNumber mate rialGroupN umber 

.I807296E+08 -.3453916E+06 -.5935107E+08 There is one line of data for each node in each element. 

.1881153E+08 .237935 1E+07 -.597 1771 E+08 


Scaling Options: 

A scaling option must be given when the data is read in so that the plotting program will know 
how to draw contours. The above data uses the fixed option which allows the user to specify the 
minim um value (-4e7 in the above data) and maximum values (4e7) when the data is read in. It is 
also possible to specify that the program automatically pick the minimum and maximum values 
when reading in the data. There are several options for doing this. These are auto, group, and 
active. 

- auto tells the program to automatically pick the min. and max. from all of the input data. 
There is no extra data necessary for this command. 

- group allows the user to specify that min. and max. value be picked from a specific 
material group. On the next line, the material group number to scale must be specified. 

- active allows the program to pick the min. and max. value from all the active elements. 
No extra data is required for this option. 

It is also possible to scale the data after it is read in by changing the Data Range fields in the 
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bottom right hand comer of the primary panel. However, the data being read in must still have 
one of the scaling options specified in the file. 

Interface 

The plotting interface consists of several Motif “widgets” or panels. These are shown in Figures 
1 and 2. The widgets shown in Figure 1, the main and redraw widgets, appear when the program 
is started. These widgets contain a collection of buttons, toggle switches, and text fields. The 
operation of each is documented below. Figure 2 shows the menu panels. The one labeled 
“Mesh” is the main panel. The “Mesh” panel is activated by using the mouse to select the main 
widget (Figure 1) and then clicking the right mouse button. On the IBM one can also press the 
F10 special function key to bring up the Mesh panel. The others are activated through the Mesh 
panel as indicated by the lines joining the panels. The menus are self-explanatory except for the 
one labeled “Modify List of Elements to be Plotted”. This panel permits one to remove a collec- 
tion of elements or to add them back. There are three methods provided for identifying the partic- 
ular elements. These are described below: 

1. Modify by Volume: Select elements whose centroids lie within the specified xyz coordi- 
nate ranges. 

2. Modify by Group Number: Select elements in the specified group. 

3. Modify by Loop List: Select elements “First” to “Last” with an “Increment” or stride. For 

example, if First, Last, and Increment are 1,10,2, respectively, then the selected elements will be 
1, 3,5,7 ,9. 


Description of Buttons, Toggles, and Text Fields on Primary Panel: 


Redraw 
Zoom In 
Zoom Out 
Node Numbers 
Element Numbers 
Shade Elements 
Label Intensity 
FontScale 
World Coordinates 


Redraw mesh using current settings. 

Zoom in on center portion of plot (magnification = 4x). 
Zoom out (reduction = 4x). 

Label nodes. 

Label elements. 

Color element according to the specified color group. 
Label element according to the specified material group. 
Magnification factor for default font size. 

Range of world coordinates in plotting window. 
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Rotation 


Magnification 
Use Displacements 


Auto World 


Rotation about z,x, and y axes - in that order followed by incremental 
rotaion about the z axis. When all angles =0, the z-axis points to the 
top of the window and the y-axis points to the right side. A right- 
handed coordinate system is used. 

Magnification factor to apply to the nodal displacements. 

Click on to plot deformed mesh. Displacements are read in using the 
menu option Displacements under Document. 

Allows program to automatically specify world coordinates for 
window based on size of mesh. 


To PostScript File Click on to create PostScript file rather than draw to screen. This 

function creates a much smaller file than saving with the default print 
command. Greyscale is always output By changing one parameter 
in file, it can be converted to color. (Directions are included in the 
PostScript file.) Currently, the output file name is “hardwired” to be 
/tmp/out.ps 


Monochrome/Color Toggles display between color and greyscale. 


Contouring 

Draw Contours 
Label Contours 

Lines 

Bands 

Outline Elements 


(All options take effect on next Redraw.) 

Click on to draw contours. 

Click on to label contour lines if just Lines selected or draw legend if 
Bands are selected. 

Click on to draw contour lines. 

Click on to draw contour bands. 

Click on to draw element bondaries when contouring. Element 
boundaries are always drawn when contouring is turned off. 


Data Range 

Min Lower limit for contour data. 

Max Upper limit for contour data. 


Comments: 

Below the Zoom Out button is an unlabeled text field. It is provided to make it convenient 
to specify a particular zoom level. For example, if the initial (ie. when the mesh file is read) hori- 
zontal and vertical ranges are 100 and 200 respectively, specifying a zoom factor of .5 will reduce 
the ranges to 50 and 100. This will result in a magnified display. 
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Chapter: 5 

Installation of Source Code and Samples 

The distribution media contains the following four compressed tar files: 
pwmeshge.z 
flex94.z 
plot95.z 
samples .z 

After copying these files to a UNIX computer, these files must be renamed as 

pwmeshge.Z 
flex94.Z 
plot95.Z 
samples. Z 

The files can then be uncompressed using the command 
uncompress * 

Next the files in each tar file are extracted using the commands 

tar -xvf pwmeshge 
tar -xvf flex94 
tar -xvf plot95 
tar -xvf samples 

The following four sub-directories are created in the carrent directory: 

PWMeshGen 

Flex94 

plot95 

Samples 


Creation of Executables 

Change to the directory containing the four sub-directories listed above, then execute the 
following commands. The words in italics are comments, not commands. 
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Flex94 


cd Flex94/Control 
make 

The executable is named “fe” and is located in the current directory. It may be moved to 
any location desired. 
cd .. 

PWMeshGen 

cd PWMeshGen 

C han ge the pathname for PWMESHGEN _DIR in the file makefile. def to 
the pathname for present working directory. 

For example : PWMESHGEN _DIR = /home4jdw/Scratch/sriren4/PWMeshGen 
Type the command: pwd to obtain the pathname for 
the present working directory. 

Similarly, change the pathname for PWMESHGEN in the file PWForm to 
the pathname for present working directory. 

For example : PWMESHGEN = /home4jdw/Scratch/sriren4// > WMer/iGe7: 

make 

The executable is named “ PWMeshGen ” and is located in the current directory. 

The executable should not be moved. 
cd .. 

Plot95 


cd Plot95 

Change the pathname for PLOTTER_DIR in the file makefile.def to 
the pathname for present working directory. 

For e xamp le : PLOTTER_DIR = /home4jdw/Scratch/sriren4/Plot95 
Type the command: pwd to obtain the pathname for the present workinz directory. 

make all 

The executable is named “ plot95 ” and is located in the Plot95 director v 


Sample Input 


Input and output files for six problems are included. These are in the subdirectories Sam- 
ple 1 - Sample6. Comments are included in the subdirectories which describe earn sample. 
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